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1.0 INTRODUCTION 


Aircraft have been designed with the NASA- developed supercritical wing to 
fly at cruise Mach numbers approaching one. The need for low-installed 
drag and high-drag-divergence Mach number nacelle installations is necessary 
to the success of this design. Design techniques are required to evaluate 
these nacelles on an isolated basis and then on an installed or integrated 
basis. 

For this reason NASA has a program underway to provide design information 
for low-drag, high-drag-divergence Mach number isolated nacelles suitable 
for use with advanced high bypass ratio turbofan engines. One element of 
such a program was the development of a method to predict the inviscid pres- » 
sure distribution and flow field about an arbitrary axisymmetric ducted body 
at transonic speeds. The prediction technique will provide the means to con- 
duct parametric studies so that the nacelle design criteria could be evaluated 
to select configurations for further experimental investigations. The predic- 
tion technique would provide guidance during wind tunnel testing to develop 
nacelle shapes which would minimize drag within given design restraints. 

Several techniques of solving the inviscid equations of motion about arbitrary 
two- or three-dimensional bodies at transonic speeds are presently available; 
however, there are no computer programs available which treat air inlet or 
nacelle configurations. The objective of the development of this computer 
analysis was the prediction of flow fields about isolated nacelles at transonic 
conditions. The solution technique was further specified to give accurate re- 
sults consistent with the requirement of relatively short computing time per 
input case as compared to that required for a time dependent finite difference 
method of solution. The method utilized to compute the flow field is the 
Streamtube Curvature Relaxation technique. 

The Streamtube Curvature Method (STC) of solving external flows has not been 
discussed significantly in the literature; however, the method is a very 
natural one. For example, engineers frequently rely on one- dimensional com- 
pressible flow relationships for a first-order solution to ducted flows. 
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The STC approach is similar except that a number of confluent streamtubes, 
with slightly different properties, are added together to obtain the total 
flow in the channel. Each streamtube is handled in much the same way as 
is the one streamtube in the one-dimensional problem. In the limit, as the 
size of the individual streamtubes approaches zero, the STC method satisfies 

the inviscid equations of motion exactly. 

> 

This report describes the structure of the computer program, Streamtube 
Curvature Analysis (STC), and the program usage and operation. Two types 
of information are included in this User's Manual, l) user oriented input 
sheets and output definitions along with one example case, and 2) programmer 
oriented descriptions showing program structure, program nomenclature, program 
messages and error codes, and operating instructions. The program listing is 
included as a separate document. 
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3.0 .PltOHLKM 0ESCKIPTI0N 

The $treamtul>c Curvature Analysis was formulated as a computer program to 
solve the inviscid equations of motion over a two-dimensional body (plane 
or axi symmetric) at transonic speeds. The computer program determines the 
flow field properties, streamlines, and pressure distribution over typical 
isolated nacelles and calculates the external pressure drag and the additive 
drag* (The additive drag is defined as the integral of pressure multiplied 
by the axial projection of the area taken along the entering streamtube from 
the undisturbed free stream conditions to the stagnation line on the cowl 
lip.) Ihe solution is the direct type in that nacelle shape, mass flows, 
and flight Mach number shall be the prime input data. 

The computer program is capable of analyzing the following geometries: 

a. Two-dimensional inviscid inlet problems without side spillage at 
zero or finite angle of attack. 

b. Axisymmetric zero-angle-of-attack isolated nacelle problem with: 

1. Short cowl nacelle in which fan duct air of a turbofan engine 
exhausts upstream of the primary air nozzle, and the pressure 

aV 

distribution on the aft nacelle (waist cowl) must be determined. 

2. Long duct nacelle in which the exhausts of both streams are 
confluent at the exit or mixed upstream at the exit. 

c. In all cases, the flow field may be calculated in the presence of 
a centerbody (or ramp) whose leading edge may be positioned either 
forward or aft of the cowl lip plane. 

The flow field boundaries can be located as far upstream, and laterally 
displaced from the nacelle as far as practicable, to ensure minimum disturbances 
at the boundary and as far downstream as necessary to ensure correct nacelle 
trailing edge flow conditions. The program is capable of computing the inlet 
internal flow up to the assumed location of the engine face. Ihe exhaust 
nozzle weight flow and aerothermodynamic properties of the exhaust flow shall 
be input quantities. The velocity of the exhaust shall be sonic or greater. 
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The computer analysis is capable of predicting the approximate location of 
any embedded normal shock in the transonic flow field over the nacelle. It 
was also desired that the analysis include free-stream Mach number's greater 
than one by using the approximate method of Moeckel for the bow shock shape, 
but the STC approach was not compatible with the approximate bow wave definition. 

The computer program was structured so that two separate computer source decks 
can be set up. One deck was limited by the requirement that the central memory 
storage could not exceed 70,000g (octal) locations on a CDC 6 400 and/or 6600 
computer and was thus restricted to a relatively coarse computational mesh 
network. The second deck utilized a finer computational mesh and included re- 
finements in the stagnation line region and in regions of high velocity grad- 
ients. The refined version was limited to a central memory storage not to ex- 
ceed l4o,OOOg (octal) locations. In actuality, the most useful deck utilizes 
768 grid intersection points and uses a central memory storage of 106,000 
octal locations. 

The size of central memory storage required is a function of the defined table 
sizes. These table sizes may be changed to meet the user’s needs. The method 
to change the table size is included in this User's Manual. Otherwise, the 
basic logic of the STC program is identical and all the capabilities defined 
above are always available. 

The coding of the computer program meets the following requirements: 

1. The program has been written to run on the LRC’s CDC GkOO and/or 
6600 computers, or any similar CDC 6400 and/or 6600. 

2. The bulk of the computer program has been coded in CDC FORTRAN 
2.3 language. At least one subroutine has been coded in CDC 
Compass 1.1 language. Ihe programs have been written to run 
under the SCOPE 3-1 operating system. 

3- All input/output has been accomplished with CDC FORTRAN 2.3 state- 
ments. The standard system file names of INPUT for card reading 
and OUTPUT for printing have been used. In addition, input from 
tape files and output to tape files has been used. 
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lhe description of the required capabilities of STC do not include all possible 
features that are included in the computer analysis. The user has control over 
the amount of grid refinement or computational mesh size, both by specifying 
local areas of mesh refinement and by setting the number of overall flow field 
refinements. The input geometry may. be specified as coordinates only or co- 
ordinates and local surface angles on any boundary. Multiple channels may be 
defined so that inlets with acoustic splitters can be analyzed. 




8 



4.0 METHOD OF SOLUTION 


4.1 Uasic Equations 

The STC Program is designed to solve the equations of motion along stream- 
lines, 'll = constant lines, and along lines which are orthogonal to the 
streamlines, 5 = constant lines. The variable § is introduced to avoid con- 
fusion with the velocity potential & which is only applicable when the flow 
is irrotational. 

Across the streamlines, the continuity and Crocco form of the momentum 
equation are written: 

Q]\l 

Continui-ty: r)A = ^ ( ? = Const) (l) 

Momentum: 1 d(V ) _ _V_ <9H _ j <dS ( ? = Const) (2) 

2 rm + oh dn 


Along the streamlines the following forms of the energy and momentum 
equations apply: 


Momentum : 



( i|> = Const) (3) 


Energy: 


DH 

Ds ~ 


0 


( ^ = Const) (4) 


where : 


A = Flow cross-sectional area = f 'JFrQn 

C = Curvature of the streamline 

H = Stagnation enthalpy 

n = Distance along the orthogonal 

p = Static pressure 

r = Radial coordinate 

s = Distance along the streamline 

S = Entropy 

T = Static temperature 

V = Velocity 

tjl, W = Stream function, cumulative flow rate 
P = Density 

The solution method is an extension of the streamline curvature method. It 
may be briefly described as follows: First a crude grid of streamlines and 

orthogonal lines are assumed. (Please refer to Fig. 4-1). Second, the 
curvature of the streamlines at each of the grid points is evaluated. Third, 
the momentum equation is integrated along a line normal to the streamlines to 
obtain velocity and the continuity equation is integrated to determine the 
"correct" streamline positions (for the assumed curvature field). These are 
indicated by the "x" in Fig. 4-1. Fourth, an adjustment, £>n, is computed by ; 
considering (l) the difference between the computed and assumed streamline 
positions and (2) the effect of the implied curvature modification in the in- 
tegrated momentum equation. Finally, the streamlines are repositioned by the 
6n values. 

Because the movement of any one grid point alters, through a change in 
curvature, the velocity at nearby points, it is highly desirable to account 
for these interrelating point adjustments simultaneously. Ihe utilization 
of a simultaneous solution procedure, employed here, is not part of the 
classical streamline curvature method [l, 2, 3^. In comparison, the classical 
method yields calculation times which are very slow, especially for a closely 
spaced calculation grid. In concept, the set of simultaneous equations for 
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2. Evaluate Curvature 


3. Integrate The Cross- 
stream Momentum Eqn. 
and Continuity Eqn. To 
Determine The "Correct" 
Streamline Positions. 


4. Solve The Matrix 

Equation for 6n and 
Move The Grid Points. 



Fig. 4-1 Solution Technique 
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the normal streamline adjustments are formulated from the finite difference 
equivalent of the following equation: 


d J (fen) + (1 - M 2 ) 
<M l 2 (PV) 2 


, = F 

as 2 


(5) 


where: 

6n = Required streamline adjustment in the normal direction 
'l* = Stream function 

s = Curvilinear distance along a given streamline 
M = Mach number 
PV = Flow per unit area 

F = Driving (or error) function derived from the solution to the 
integral continuity and normal momentum equations. 

This equation is derived in (Ref. 4) for the special case of isentropic 
2-dimensional flow. (These limiting assumptions are utilized only to main- 
tain simplicity of illustration; they are not part of the computer program. ) 
From a mathematical point of view, the above equation is similar to the con- 
ventional equations for velocity potential or stream function, namely, 


(1 - 


M 2 ) 

y 


5*~0 

ay 2 


M 

y ay 


♦ ( 1 - m x 2) £-fL. 0 

ax 2 


(6) 


(i - m 2 ) 5 2 * 2M m a 2 '!' (i - m 2 ) a 2 '!' 

y — - - x y - — - — + x “r = o 


ar 


9x ay 


(7) 


0x 


where 



x and y are coordinates in a two-dimensional rectangular system and 0 and 
are the velocity potential and stream function, respectively. (Again, all 
of the equations have been restricted to isentropic flow for illustration. ) 

However, for the purpose of calculating transonic flow, the use of Eq. (5) 
offers a distinct advantage over Eqs. (6) or ( 7 ) for the following reason: 
Because the grid is always alligned in the streamline and normal to stream- 
line directions, no cross-derivatives appear. In consequence, the finite 
difference star is simply switched from a subsonic representation to a super- 
sonic representation illustrated in Fig. 4-2. 

Notice that for supersonic flow, no points downstream of the orthogonal line 
are included. This reflects the physical reality that disturbances downstream 
will not affect the flow upstream. It is, of course, because the coefficient 
term, (l-M 2 ), passes through zero and changes sign that the star-switching 
noted in Fig. 4-2 is appropriate. 

If the grid system is not alligned with the flow direction, a cross derivative 

appears in the equation, as in Eq. (6). Unfortunately, the mixed derivative 

coefficient, M M , does not have the same sign change property and, therefore, 
x y 

star-switching with this equation is more difficult for general flows with 
some angularity. 

The extended streamline curvature method, here referred to as the Streamtube 
Curvature (STC) method, has then the following features: 

o No additional complexities arise when the flow is rotational. 

o The slip lines between the exhaust jet and the external flows can 
be handled precisely. (The procedure is to consider two coincident 
streamlines. Their position and pressure are the same; their ve- 
locity and stagnation properties may be different.) 

o From a numerical point of view, the streamline/orthogonal line oriented 
grid facilitates the analysis of transonic fields, as described above. 

o The streamline/orthogonal line grid also provides a mapping of the 

flow field into a rectangular domain. This is helpful from the stand- 
point of computer program organization. 
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The STC Program has also been designed 
o to handle multiple streams, 

o to adjust the flow rates of the jet exhaust streams to meet the 
(2-dimensional) choking condition, 

o to place grid points at locations in the flow field where they are 
needed, as determined by local variations of the dependent variables, 
and 

o to allow external flow analysis by incorporating matched near- field 
and far-field solutions. The far-field solutions are obtained 
analytically, utilizing small perturbation theory. 

/' 

4.2 An Outline of The Calculation Steps 

The operations performed by the STC Program may be outlined as follows: 

1. Define tile flow regions and locate (approximately) the "primary" orthog- 
onals and the streamlines which divide the internal and external flows. 

2. Refine the grid as required by inserting additional streamlines and 
orthogonal lines between those already existing. 

i. 4 * ' 

3- Compute the streamline angles and curvatures. 

4. Compute the orthogonal line angles and move the grid points along 
the streamlines to obtain orthogonality. 

5. Compute the velocities on the "far-field" boundary. 

6. Adjust the flow rates in the exhaust streams, if any, to meet the 
calculated choking flow rate. 

7. Integrate along each orthogonal the momentum and continuity equations 
C Eqs. (l) and (2) 3 . 

8. Determine if the streamline positions are within a "rough tolerance". 

If so, return to Step 2 for additional grid refinement (unless grid re- 
finement limits have already been reached). Otherwise, continue to 
Step 9* 
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'). Do tormirio if the .streamline* positions are within final tolerance. 

If mo, jump to Step 13. Otherwise continue to Step 10. 

10. Set up the matrix equation for the streamline correction, &n. 

11. Solve the matrix equation. 

12. Modify the streamline positions by 6n, and return to Step 3. 

13. Calculate and print the output quantities; then return to Step 1 
for the next case, if any. 

The first operation includes reading the card input for a description of the 
geometry and flow properties. The computer program has been written to have 
general capability for analyzing a great variety of configurations. The 
first step in the programmed logic is to develop a table of orthogonals or. 
calculation stations for each of the several flow regions. The regions are 
determined as illustrated in Fig. 4-3 so the calculation can proceed from up- 
stream to downstream. The boundary of each region is defined as a primary 
orthogonal. As shown in Fig. 4-4, the initial grid which is developed con- 
tains only the primary orthogonals and the double streamlines which separate 
the various streams. 

The second step in the computational procedure is the grid refinement. The 
very crude grid, obtained in Step 1, is refined before the first solution of 
the flow field equations is executed. A new orthogonal is placed within each 
region and, likewise, a streamline is inserted in the middle of each channel. 
In the external channel, additional streamlines are placed close to the body. 
After the solution has been obtained for this net, the grid intervals are 
halved as required. This may be likened to the steps taken when one "flux 
plots" a flow field by hand. First, major flow lines and normals are sketched 
in, and then more and more streamlines and orthogonal lines are added until 
the desired resolution is obtained. At each step in the process the positions 
of the lines are adjusted to meet the correct solution requirements. Ihe pro- 
cedure automatically provides for grid refinement in regions of high curvature 
and high acceleration or deceleration. The streamline and orthogonal lines 
which are added between existing lines are not required to span the field if 


16 




17 


FIGURE ,4-3 SUBDIVISION OF THE FLOW FIELD INTO REGIONS 



FIGURE W INITIAL STREAMLINE/ ORTHOGONAL GRID BEFORE THE FIRST REFINEMENT 


i 

only local rol’i demerit, noar the body, is required. The refinement procedure 
presently' built into the program uses a criteria involving the distance and 

j?' 

„•/*' . • - 

velocity increment between grid points. These refinement criteria are dis- 
cussed in detail in Section 10 and 12. ' 

The third step in the method is to determine the angles and curvatures of 
the streamlines at each grid point. Fbr subsonic portions of the flow field 
this is performed by fitting a piecewise continuous cubic polynominal in a 
coordinate system which is locally rotated for each interval. The resulting 
fit is analogous to the curve produced by a beam which is loaded by discrete 
forces so as to pass through the given grid points. The locally rotated co- 
ordinate system removes the restriction that requires the slope to be small. 
For grid points located in a supersonic region, backward difference formulas 
are employed. Either 3-point, 4-point, or 5-point formulas may be optionally 
selected. Again the coordinate system is rotated so that slopes in the curve 
fitting coordinate system are small. 

In the fourth step the orthogonality of the grid points is checked and points 
'are moved along the spline curve as required to achieve normal intersections 
between the two sets of lines. Also, the normal distance, n, is computed for 
each grid point as measured from the lower boundary of the orthogonal. 

s. ■ 

When the initial grid is set up, a boundary is placed some distance away from 
the body. This boundary becomes the interface between the near-field and 
far-field solutions. The near-field is computed by the streamtube curvature 
method and the far-field is computed by linear small perturbation theory. In 
the process of iterating, this boundary (which is also a streamline) will 
float so that its shape and velocity distribution are matched by both the 
inner and outer solutions. In practice, the shape of the interface stream- 
line (also referred to as the far-field boundary) is first assumed. Using 
the far-field equations, the velocity distribution is calculated. This is 
Step 5. These velocities are subsequently employed in the near-field analysis 
and from this comes a revised shape for the Interface streamline. Revised 
velocities will then be computed in Step 5 during the following iteration 
cycle, and so forth. 
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Step 6 is the modification, as required, of the flow rates of the exhaust 
streams. For boattail analysis of nacelles, the internal geometry of the 
exhaust passage is required input to the SIC program. Because of streamline 
curvature effects, the discharge coefficient for the nozzle will be somewhat 
less than unity. The user, however, may input a flow rate based on unity 
discharge flow coefficient or, for the matter, any approximate value. De- 
termination of the velocity distribution across the throat of the nozzle 
will be determined within the STG framework and the evaluation of the maximum 
"choked" flow rate is Step 6 of the calculation procedure. 

Step 7 is the solution of the flow field equations per se. This section of 

the program is referred to as the "flow balance"; Eqs. (l) and (2) are inte- 

! 

grated. In the external regions of the field, the momentum equation is in- 
tegrated from the far-field interface boundary to the body (or to the center- 
line or lower boundary, whichever exists). Hie integral form of the momentum 
equation is: 

n FF 

ln V o 2 - ln V - 2 L c dn (8) 

o 

whore: 

V =. Velocity as determined in Step 5 along the far-field streamline. 

1* r 

= Velocity at any streamline with orthogonal distance n^. 

n = Distance measured along the orthogonal to the far-field 
FF 

streamline. 

Although not reflected in Eq. (8), the effect of varying total pressure 
behind a shock wave is also included. Also, if a slip-line occurs in the 
field, the velocity jump equations, 
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(9a) 





2 


H 


(9b) 



(9c) 

i 


t 


are employed where the subscripts (+) and (-) denote conditions on the 
streamlines above and below the slip-line, respectively. 


The velocity, total temperature and total pressure allow determination of 
the density at each grid point, and the inverse product of density and ve- 
locity is integrated to find flow area. 

r 



eft 

Pv 


do) 


The cumulative flow areas calculated by Eq. (10) are compared with the 
geometric areas of the streamlines used in Step 3- The difference between 
these two values is used as a convergence check (Steps 8 and 9) and in the 
streamline correction equation, Step lO. 


For internal flow orthogonals, the velocity at the outer boundary 

(V in Ea (8) ) is not known. Instead, an iteration process is employed 
FF 

whereby the outer boundary velocity is varied to obtain a match in the. cal- 
culated geometric passage area. 

In Steps 10 and 11 the proper adjustment of the streamline positions is de- 
termined and in Step 12 the grid points are moved in the normal direction 
by this computed adjustment. 
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The iterative sequence is to start with a crude grid, as noted above, and 
to repeat Step 3 through 12 until the flow balance error is small. This is 
often accomplished in one or two iterations. The grid is then refined to 
the next level and the field is reconverged. The refinement/convergence 
process is continued until the grid refinement criteria is satisfied, or 
alternately, until computer storage limits are reached. At this point, ad- 
ditional loops through Steps 3 to 12 may be performed until the flow balance 
error is satisfactory. 

A complete description of the details of the numerical procedure is contained 
in Reference A. 



5.0 PROGRAM DESCRIPTION 


5.1 Data Storage Tables 

The framework of the STC program is designed to allow flexibility as to 
the configuration to be analyzed. For example, very weak limits are placed 
on the number of flow boundaries and the number of channels into which the 
flow is split. And no specific limits' are placed on the number of stream- 
lines or the number of orthogonal lines in any given region of the flow. 

To acomplish this, the bulk of the calculation data is saved in arrays which 
are singly dimensioned. Within each array the data are packed together 
to maximize storage efficiency. Descriptions of these tables are provided 
below. 


5.1.1 Field Tables 


Flow field data are stored in singly subscripted arrays, one for each 
variable. The quantities saved at each streamline -orthogonal line grid 
intersection are as follows : 


C. 


Z 

R 

VM 

PHI1 

CURV 

51 

52 
B 


RHS 

DS2 


axial position, abscissa 

radial or vertical position, ordinate 

velocity 

streamline angle, measured from horizontal 
curvature of the streamline, -d<j>/ds^ 
curvilinear distance along the streamline 
curvilinear distance along the orthogonal line 
coefficient in matrix equation for DS2, 
indicator of subsonic (B>0) or supersonic 
(B<0) velocity 

right hand side of matrix equation for DS2 
correction of streamline positions 


In the Fortran coding M is the symbol commonly used as the subscript 
for the above arrays and NM is the number of field points (or maximum 
value of M) . The points are grouped by orthogonals starting at the up- 
stream orthogonal. Points along the orthogonal are grouped together. 
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Hence, the point below and above the point are (M-l) and (Mfl) res 
pectively. The neighbors in the streamwise direction, however, are 
determined by referring to the JMS-table described in Section 5.2. 

5.1.2 Channel Input Data Table 

Information such as the boundary coordinates and flow properties is 
compactly stored in a single array, TABLES, so that only the total 
amount of information saved is limited by the array size. No limit is 
placed on the amount of information to be placed in any one Table of 
which there are six: 

o channel input table 

o boundary table 

o table of convected properties 

o table of wake displacement thickness 

o flow adjustment table 

o station table 

In the first table, CHDATA, input information read from page STC/ 
sheet-3 is stored. The information is stored in subtables, one sub- 
table for each channel, and the arrangement of each subtable is as 
follows: 


CHNAM 

channel name (BCD) 

LHNEXT 

length of subtable 

WTFLOW 

flow rate (if input) 

TTO 

total temperature 

PTO 

total pressure 

TSO 

static temperature 

PSO 

static pressure 

MACHO 

Mach number 

AO 

flow cross-section area 

VARY 

indicator as to whether the 
flow rate may be varied 

RG 

gas instant 

GAM 

ratio of specific heats 
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NR not operational 

NC " " 

TAB " " 

BB(NR&NC) not operational, zero length array 

Except for the first two words, CHNAM and LHNEXT, all of the above input items 
are optional input and they are all equal to BITS unless values are supplied 
according to input instructions for STC/Sheet-3. If data for a second chan- 
nel are supplied, these data will follow the first channel and so forth 
for any number of channels. The first word of the CHDATA table is at location 
LHO and the last word is at location LHE relative to the origin of the 
TABLES-array . If LHE=LHO-l, no channel data has been input and the channel 
data table has zero length. In the Fortran coding, the subscript LH is 
used to refer to the channel data. 

Channel information is read into the STC program by subroutine RCD; the 
stored channel data is utilized by routines RTCFI, BCONV, ADJWF and ISBOT. 

5.1.3 Boundary Table 

Directly following the channel data are the coordinates of the boundaries. 
Again a subtable for each boundary is constructed and the information is 
stored in the following order: 


BDT 

boundary name 

LBNEXT 

length of subtable 

LBZ1 

index increment to first coordinate 
in ZBT, RBT, ANGBT-lists 

CHNAME 

channel with which the boundary data 
is associated 

UP 

upper or lower boundary indicator 
(if the boundaries around a leading 
edge have been collated together, then 
CHNAME ( 2) =UP is the name of the channel 
above the leading edge and CHNAME(l) is 
the name of the channel below the leading 
edge) . 
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LEDEX 


ZBT(l) 

RBT(l) 

ANGBT (1) 

ZBT(2) 

RBT(2) 

ANGBT(2) 
o 
o 
o 
o 

BDNAME name of a specific boundary when several 

boundaries are collated together in one 
subtable 

LBA, LBB index limits relative to ZBT, RBT, ANGBT of the 

coordinates for BDNAME 

(NOTE - the above 3 items are repeated for each boundary 
where LBZ1/3 is the number of collated boundaries. The 
existence of this information results in the displacement 
of the ZBT, RBT, ANGBT coordinates to higher memory 
locations. BDNAME, LBA and LBB are equivalenced to ZBT, 

RBT and ANGBT respectively) . 

Boundary coordinates are necessary input to the program so two or more 
subtables will always exist, stacked one after another. The first word 
of the BDYTAB table is at location LBDO and the last word is at location 
LBDE. LB is the subscript used to refer to the boundary data. 

Boundary coordinates as supplied on page STC/Sheet-2 of the input sheets 
are read by routine RBD. This routine converts the angles from degrees 


index (relative to ZBT, RBT, ANGBT) of the 
leading edge point when boundaries are 
collated 

axial coordinate (x) 

radial coordinate (y) 

surface angle (measured from x-axis in 

radians) 
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to radians, translates and rotates the coordinates as required and stores 
the points so that as one proceeds from one point to the next (or walks 
along the contour) one's left arm is next to the flow field. Thus, the 
ordering of points on an upper flow boundary will be reversed and the 
angles incremented by 180 degrees. 

The several boundaries which may be defined to comprise one continuous 
contour are collated in subroutine BLDTAB. The boundary table is 
referenced in routine BDYPTM. 

5.1.4 Table of Convected Properties 

For each channel, a subtable of convected properties and channel flow 
data is built by subroutine BCONVT from data in the channel input table 
(if it exists) and input data from STC/sheet 1. This table contains some 
of the same information as the CHDATA table. Information contained in 
the CONVTB table is complete and follows a consistent format. 


CH 

channel name 


LTNEXT 

length of subtable 


NPT 

not operational (“1) 


LPSI 

(=15) 


LTT 

(=16) 


LPT 

(=17) 


LRCU 

(=18) 


CRG 

gas constant 


CPGJ 

specific heat at constant 

pressure 

C2CP 

twice CPGJ 


QGAM 

inverse of ratio of specific heats 

FGT 

(y-D/y 


FGP 

Y/(Y-1) 


FGR 

1/(Y-1) 


AREATB 

area for calculating flow 

rate 

PS I 

flow rate for the channel 
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TT stagnation temperature 

PT stagnation pressure 

RCU not operational 

Again the data for the several channels are stacked one after another. 
LTO and LTE are the first and last locations of the table relative to 
the origin of TABLES. LT is the subscript used to retrieve convected 
property information in subroutines BLDTBS, RBCONV, ADJWF and TTPT. 

5.1.5 Table of Wake Displacement Thickness 

The wake displacement table is only constructed if there, exists a 
blunt trailing edge, and is arranged as follows: 


X2W 


streamline coordinate, 

LWNEXT 

length of subtable, = 2N+2 

S1W 

(1) 

list of distances from trailing edge 

S1W 

(2) 


S1W 

(N) 


DST 

(1) 

list of wake displacement thicknesses 

DST 

(2) 


DST 

(N) 



This table is built in the subroutine BLDTBS. DST is determined so that 
the wake thickness equals the trailing edge thickness at S1W=0. 
Thereafter the wake thickness decrease at the rate of 0.1 times the 
distance from the trailing edge. At 10 trailing edge thicknesses down- 
stream, the wake thickness is zero. 

As many wake thickness subtables are built as there are blunt trailing 
edges. The wake thickness table begins within the TABLES-array at 
location LWO and ends at LWE. The table is referenced by the TTPT 
routine. 
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5.1.6 Flow Adjustment Table 

The flow adjustment table is also created in the BLDTBS routine, one 
subtable for each trailing edge. The information contained is as 

orthogonal coordinate of the t.e., 
streamline coordinate of the t.e., £2 

coordinate of the choked station of the flow 
below the t.e. if not X1F 

coordinate of the choked station of the flow 
above the t.e. if not X1F 

curvilinear streamline distance to the trailing 
edge (along the upper surface of the airfoil) . 

This value is used for interpolating the wake 
displacement thickness, 
number of channels below the t.e. 

" " " above " " 

= -1 if the single channel flow is choked 

* 1 if channel flow rates below the t.e. are known 

= 2 if channel flow rates above the t.e. are known 

12 element storage array used by subroutine 
NEWRAP for the flow iteration 

The data stored in this table is used by subroutine ADJWF which adjusts 
the flow so that at each trailing edge the pressure difference from one 
side of the trailing edge to the other is reduced to zero (to satisfy 
the kutta condition). If this condition cannot be satisfied, the flow 
in one of the channels will be choked and routine ADJWF adjusts the flow 
to the maximum choked value. The flow adjustment table is located be- 
tween LFO and LFE. Each subtable is NFCOLS (°20) in length. 

5.1.7 Station Table 

The station table is the last of the compacted tables (except for the 
shock point table which is appended to it) and it contains information 


follows : 

XI F 
X2F 
X1BF 

X1AF 

S1F 


NCHB 

NCHA 

JORDER 


VNR (12) 
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for each orthogonal line. The table grows during the calculation process 
because the number of orthogonal9 is Increased to obtain a refined grid. 
Because it is the last table it can easily be extended into the unused 
portion of the allotted memory. The data saved for each of the orthogonals 
is arranged in the station table as follows: 


XI 

LNEXT 

MLB 

MUB 

PRIM 

TYPELB 


NAMELB 

1LB 

FLB 

S1LB 

TYPEUB 

NAMEUB 

IUB 

FUB 

S1UB 

VMB 


station coordinate, 
length of subtable 

field point index of first (lower boundary) 
point of the orthogonal 

field point index of the last (upper boundary) 
point of the orthogonal 

primary station indicator, T or F, (a primary 
station, is one of the original grid stations 
and it will pass through boundary end points), 
type of the lower boundary, i.e. SOLID, 

FIELD, FARFLD, etc. for indicating proper 

boundary condition 

name of the lower boundary used for 

referencing the boundary table 

boundary (table) interval of the orthogonal 

boundary intersection 

fractional position in that interval 

curvilinear distance from beginning of the 

interval 

type of the upper boundary 
name of the upper boundary 
interval of the upper boundary 
fractional position in the interval 
curvilinear distance from the beginning 
of the interval 

boundary velocity, used as an initial guess 
for the velocity iteration in the FLOBAL 
routine 
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DWDV .SCHOKE 


X2CL 


SLSWI 

MCL 


slope of the flows rate versus velocity 
curve used for the velocity iteration in 
the FLOBAL routine, = SCHOKE if flow is 
choked at this station 

coordinate of the control streamline 
used for positioning orthogonal lines 
in the PTMOVE routine 

sonic line-shock wave indicator, =1.0 for 
mixed flow, = 0. otherwise 
field point index of control streamline, 
used in PTMOVE only 


This group of 20 items is repeated for each station, starting with the 
upstream and proceeding to the downstream stations. Field point 

^calculations are then performed by looping through the station table, 

/ 

starting with the first station at L=L0 and proceeding to the last 
station. The last word of the last station is at location LESTA. 

’The above set of six tables is initially constructed with gaps between 
each table. The channel data and boundary tables are constructed simul- 
taneously in the RCD and RDC routines called from REDINP. The maximum 
size of the channel data routine is preset by the value of MAXLH. 

MAXLH is initialized to 400 but may be input as any value if such should 
be necessary. 


After reading the card data input, at the beginning of the BLDTAB routine, 
the boundary data is moved down in memory so that it is just above the 
channel data. The gap between the two tables is thus eliminated. 

The remaining four tables are then all built in subroutine BLDTBS . Again 
spaces between the tables are provided so that they can be constructed 
simultaneously and grow to the following lengths without interfering with 
each other: 
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Table Name 


Input Variable Default 

Name Value 


Convected Properties M/UCLT 200 

Wake Displacement Thickness MAXLW 200 

Flow Adjustment MAXLF 200 

At the end of the BLDTBS routine the gaps between the tables are eliminated 
by moving each table down in memory. The result is that all information 
contained in these tables is stored compactly, and the last table, the 
station table, has space into which it may expand. The length of all six 
tables combined is limited to the length of the TABLES-array , and this value 
is adjusted at program load time to meet the size requirements of the flow 
field. 


5.1.8 Streamline Table 

The streamline table is not stored in a compact arrangement. It consists of 
three arrays in a la bled common SLTAB. 

W cumulative flow rate for the channel 

X2 % - coordinate 

SLCHN channel name of which this streamline is a part 

The streamline number, J, is the subscript used to access information in 
each of these three arrays. All of the streamlines for a given channel are 
together and in order (proceeding away from the centerline). However, no 
special ordering of channels is required. 


5.1.9 Table of Leading Edge and Trailing Edge Points 

The Leading Edge- Trailing Edge Point (LETEPT) table is constructed in routine 
BLDTAB and used in routine BLDTBS for the purpose of defining the flow regions 
and primary orthogonals. The data in this table is not saved after the BLDTBS 
routine is executed. 

The information in this table is obtained from the boundary table. The ends 
of boundaries and double points contained within the boundaries are listed, 
together with the boundary and channel names. Each "line" of the LETEPT table 
contains the following ten items. 


32 



axial coordinate of boundary end point or 
double point 

radial or vertical coordinate 
mean angle of the flow at XE, YE 

number of times the same point has occurred as an up- 
stream boundary end point. Normally NLE=0 or 1. If 
NLE=2 then the point is a leading edge, 
number of times point has occurred as a downstream 
end point. NTE=2 for trailing edge point. NLE=NTE=0 
for a double point in the boundary table 
name of flow stream above point 
name of flow stream below point 1 
name of boundary (UPPER=F) 
name of boundary (Upper=T) 

number of times the point has been used in developing 
the ORTCHN- table, initially = -1 for a double point 

After the table is constructed the points (or lines) are sorted so that the 
upstream points are first and the points follow from upstream to downstream 
in order. If difficulty is encountered in the development of the basic grid, 
the information contained in the LETEPT-table may be helpful in diagnosing 
the error. 

5.1.10 Table of Channels Embraced by Each Orthogonal 

The LETEPT-table contains all points through which (primary) orthogonal lines 
are to pass. From this table the ORTCHN-table is developed as an aid to the 
construction of the initial grid of streamlines and orthogonal lines. The 
latter table contains a list of all of the channels embraced by each orthogonal. 
Specifically each "line" of the table contains: 

LEDGE index of the point in the LETEPT-table 

LRPREV previous line number (or ORTCHN-table index) 

of the upstream orthogonal 
CHNA(l) channel names 

CHNA(2) " " 

CHNA » " 

CHNA(2) " " 


XE 

YE 

ANGE 

NLE 


NTE 


CHL 

CHU 

BDL 

BDU 

NUSED 
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The above information is tabulated (one line) for each orthogonal. Also, 
the first two lines are developed as dummies for the purpose of listing all 
channels. Other dummy lines may also be inserted. LRD is the index in- 
crement between lines and is equal to the total number of channels minus two. 

5.1.11 Table of Index Limits 

The labeled common /IX0RIG/ contains the index limits for each of the above 
listed tables. The items and order of storage are as follows; 

LHO location of channel table origin 

LHE location of channel table end 

LBDO location of boundary table origin 

LBDE location of boundary table end 

LTO 1 location of convected properties table origin 

LTE location of convected properties table end 

LWO location of wake displacement thickness table origin 

LWE location of wake displacement thickness table end 

LFO location of flow adjustment table origin 

LFE location of flow adjustment table end 

LO location of the station table origin 

LESTA location of the station table end 

LSO location of the shock point table origin 

IJ5E location of the shock point table end 

LDUM(6) unused 

MO unused 

NM number of field points 

NJ number of streamlines 

NFCOLS number of "columns" in the flow adjustment table 
MAXNJ maximum number of streamlines (dimensional limits are MAXNJ-128) 

MAXOL maximum number of points on any one orthogonal (dimensional 
limits are MAXOL=96) 

MAXNM maximum length of field arrays, calculated in subroutine REDINP 
MAXLE maximum table length, calculated in subroutine REDINP 

LEO location of the first word in the LETEPT-table, (=l) 
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LEE 


LRO 

LRE 


location of the last word in the LETEPT- table 
location of the first word in the ORTCHN-table 
location of the last word in the ORTCHN-table 
the ORTCHN-table is subdivided into lines and LRD is the 
length of the lines 


LRD 



5.2 Accessing Data in the Field Tables 

The data stored in the field tables consist of the Z, R, SI, S2, PHIl, 
CURV, VM, B, RHS, and DS2 at each point in the field. An additional 
array termed the JMS table provides access to the information in the field 


tables. This 

information (subscripted M=1,NM) consists 

of : 


J 

' • <; 

Streamline number 



MU 

M subscript of upstream grid point; MU = 
of streamline. 

0 at beginning 


MD 

M subscript of downstream grid point; MD 
streamline . 

* 0 at end of 


ISTAG 

Point type indicator 
ISTAG = 0 - Normal point. 




1 - Stagnation (or singularity) point. 

2 - Trailing edge point or a point fixed on the 

body surface used to locate a primary, 
orthogonal . 

3 - Point adjacent to a stagnation (singularity) 

point or an end point of a partial orthogonal. 


v ‘ The word 

content of each JMS 

entry is: 




- 

GE-635 







(36 bit word) 


J 

MU 

MD 

ISTAG 




8 

13 

13 

2 




bits 

bits 

bits 

bits 


CDC-6400/6600 







(60 bit word) 


J 

MU 

MD 

ISTAG 



24 

8 

13 

13 

2 



bits unused 

bits 

bits 

bits 

bits 



On the GE-635, the machine word length presents a practical limit to both the . | 

i 

number of streamlines (J = 255) and the maximum field size (NM - 8191) It j 

was felt that these limits represent a convenient maximum for any problem ! 

which might be encountered. Hence, the bit configuration for the JMS word j 

on the CDC 6400/6600 is identical to that on the GE635, resulting in an 
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unused por.ti.on of 24 bits. The subroutines which pack and unpack the JMS 
table entries (SAVIX, GETIX, and GETRLX) are coded in the CDC assembly 
language COMPASS 1.1 in the interest of increased computation speed. These 
routines would have to be modified if a problem required a streamline number 
J in excess of 255 or a grid point index M greater than 8191* 

5.3 STC Calculation Steps. Flow Chart, Overlay Description 

The present section outlines the sequence of calculation steps performed 
by the STC Program and describes the overlay structure on the CDC 6400/6600 
computers. 

5.3»1 STC Calculation Steps & Flow Chart 

The processing flow chart for the STC program is shown in Figure 5-1 and in- 
cludes the principal subroutines, their function, and their associated output. 
The general calculation steps performed by the STC Program were outlined in 
Section 4, Method of Solution. Here the calculation steps will be identified 
with the subroutine performing that step. 

The computer program has been written to analyze many types of geometries 
and flows. The user must identify the various flow boundaries and channels 
when the program input is compiled. These sets of flow boundaries (BDY) and 
channel names (CHN) define flow regions with specified properties of temper- 
ature, pressure, and velocity. It is the task of the program user to organize 
his problem so that each flow region is defined by the proper boundary and 
channel names (see Figure 4-3 and 4-4). 

The first operation in the computer program consists of reading the card in- 
put for the description of geometry and flow properties in each flow region 
(REDINP), and then storing and building the various tables (BLDTAB, BLDTBS, 
BCONV). The channel data, boundary data, and flow property data are each 
stored in their respective tables. If necessary, the boundary points are 
smoothed and local angles at each boundary point are calculated. 
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The field point table and station table are started with the first unrefined 
grid of orthogonals and dividing streamlines. The boundary of each region 
is defined as a primary orthogonal. The dividing streamlines which separate 
the various streams are called double streamlines. 

The next step in the computational procedure is to refine the very crude 
initial grid (REFINE) after calculating the distances between points along 
streamlines. A new orthogonal is placed within each region and a new stream- 
line is inserted in the middle of each channel. In the external channel 
(identified by CHN EXT), additional streamlines are placed close to the 
body. After the solution for this grid has been obtained, the intervals are 
again halved as required. 

As the calculation procedure continues, grid refinement is automatically pro- 
vided in specified flow regions or in regions of high curvature and high ac- 
celeration and deceleration. The streamline and orthogonal lines which are 
added between existing lines are not required to span the field if only 
local refinement near a boundary is needed. These refinement criteria are 
discussed in Section 10 and 12. 

After refining the grid, the next step (SLC) in the solution is to determine 
the angles and curvatures of the streamlines at each grid point. For subsonic 
portions of the flow field, this is performed by fitting a piecewise continuous 
cubic polynomial (beam) in a coordinate system which is locally rotated for 
each interval (Figure 5-3 )• The locally rotated coordinate system removes 
the restriction that requires the slope to be small. The matching conditions 
are the angles and curvatures at each point. At the end of streamlines which 
terminate at a flow exit boundary or extend to a flow inlet, the end curvature 
is specified. Normally the end curvature is zero, but the user may input a 
constant non-zero value of curvature. 

For grid points located in a supersonic region, the subroutine SIX employs 
backward difference formulas in keeping with the switch to the star with no 
downstream points (Figure 4-2). Either a 3-point parabola, a 4-point piece- 
wise cubic (beam), or a 5-point formula may be optionally selected. The 3-point 
parabola is preferred. Again the coordinate system is rotated so that slopes 
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in the curve fitting coordinate system are small. The end conditions for 
the supersonic curve fit formulas are a specified angle and zero curvature. 

One additional task performed by SLC is the location of stagnation points 
and the definition of the dividing streamline intersecting the boundary at 
the stagnation point. At leading edges, the dividing streamline is set per- 
pendicular to the boundary surface at the stagnation point, (Figure 5-^)* 

An orthogonal which goes with the stagnation point is defined and the point 
on the first streamline is positioned. 

PTMOVE checks the orthogonality of the grid points and moves the points along 
the spline curve as required to achieve normal intersections, (Figure 5-5), 
between the two sets of lines. Also, the normal distance, n, is computed 
for each grid point as measured from the lower boundary of the orthogonal. 

When the near-field grid is defined, a boundary is placed some distance away 
from the body, lhis boundary becomes the interface between the near- field 
and far-field solutions. The near-field is computed by the streamtube curvature 
method and the far-field is computed by linear small perturbation theory, 

(Figure 5-6). In the process of iterating, this boundary (which is also a 
streamline) will float so that its shape and velocity distribution are matched 
by both the inner and outer solutions. In practice, the shape of the inter- 
face streamline (also referred to as the far-field boundary) is first assumed. 
Using the far-field equations, the velocity distribution is calculated. These 
velocities computed by FARFLD, are subsequently employed in the near field 
analysis and from this comes a revised shape for the interface streamline. 
Revised velocities will then be computed in FARFLD during the following iter- 
ation cycle, and so forth. FARFLD is a selected boundary option when the far- 
field boundary, BDY, is called FF. A solid boundary may also be specified 
by renaming the far-field boundary and defining the coordinates. 

The next step (ADJWF) is the modification, as required, of the flow rates of 
the exhaust streams. For boattail analysis of nacelles, the internal geometry 
of the exhaust passage is required input to the STC program. Because of 
streamline curvature effects, the discharge coefficient for the nozzle will 
be somewhat less than unity. The user, however, may input a flow rate based 
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on unity discharge flow coefficient or, for the matter, any approximate value. 
Determination of the velocity distribution across the throat of the nozzle 
will be determined within the STC framework and the evaluation of the maximum 
"choked" flow rate is ADJWF of the calculation procedure. 

The next three subroutines, S1AL00, BRHS, and FLOBAL, develop the solution to 
the flow field equations represented by continuity, Eq. 1, and radial or cross- 
stream momentum, Eq. 2, (Section 4 ). At each station ( STALOO) along a boundary, 
the flow equations are integrated along an orthogonal (FLOBAL) and the right hand 
side and the coefficient B (B = [(l - in Eq. 5 are determined (BRHS). 

In the external regions of the field, the momentum equation is integrated from 
the far-field interface boundary to the body (or to the centerline or lower 
boundary, whichever exists). The cumulative streamtube flow areas are then 
calculated by integrating the continuity equation and compared with the geo- 
metric areas of the streamlines defined by SLC. The differences between the 
streamline position determined by the "flow balance" FLOBAL and by SLC is used 
as a convergence check. It also defines the right hand side or driving error 
function in the streamline correction equation, Eq. 5. 

If the streamline positions are within a convergence tolerance the program 
logic loops back to REFINE. Additional streamlines and orthogonals are added 
as required by the refinement criteria. The refinement procedure presently 
built into the program uses a criteria involving the distance and velocity in- 
crement between grid points. 

If the streamline convergence tolerance has not been met, the streamline cor- 
rection equation is solved by defining a matrix, MCOEF, and by solving this 
matrix using one of several methods. The solution techniques include IAD, 
Implicit Alternating Direction, or LRELAX, Streamline Block Relaxation. The 
user has control of the choice of solution method, but IAD has demonstrated 
the more stable and faster solutions. 

After calculating the streamline correction, 6n, the streamlines are adjusted 
(ADJSL) to their new positions and the program logic loops back to SLC to de- 
fine new curvatures and geometric streamtube areas. This loop is called an 
inner solution. Depending on the amount of grid refinement, the inner loop 
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may take from 2 to 20 passes to reach streamline convergence. 

As soon above, the iterative sequence is to start with a crude grid and to go 
through subroutines SLC to ADJSL until the flow balance error is small. Ihe 
grid is then refined to the next level and the field is reconverged. The re- 
finement/convergence process is continued until the grid refinement criteria 
is satisfied, or alternately, until computer storage limits are reached. At 
this point, additional inner loops may be performed until the flow balance 
error is satisfactory. 

Finally, the output quantities are calculated and the results are printed in 
any of several forms. The output forms are defined in Section 7* 



5.3.2 STC Overlay Description ' 

The STC program has been structured for execution on the CDC 6400/6600 
machines under the SCOPE 3.1 operating system. The basic OVERLAY features 
of SCOPE 3.1 have been utilized to reduce the memory requirements to those 
currently in use at the NASA Langley Research Center. Shown in Figure 5-2 
is the overlay structure including all subroutines and the important data 
table storage areas. As indicated in this figure, the program consists 
of a main overlay, four (4) primary overlays, and six (6) subordinate 
secondary overlays. A brief description of the processing in each overlay 
is given in the following section. In this description, the word "link" is 
used interchangeably with the word "overlay". 

Overlay (0, 0) - Entry STCA 

The main overlay (0, 0) contains the main program STCA, as well as 
general purpose subroutines which are called by subprograms in the subse- 
quent primary and secondary overlays . The main program STCA provides 
control for execution of the STC program by loading the appropriate over- 
lays in sequence and testing for completion or convergence of the problem 
solution. 

Also included in the main overlay is the data storage area consisting 
of the general program tables, the field tables, and the streamline table 
as discussed in Section 5.1. This storage area represents the largest 
block of working storage in the main overlay. The sizes of these various 
regions of memory are initialized by the block data subprogram USECDG. An 
increase in the number of available field points, stations or streamlines 
may be effected simply by recompilation of this block data subprogram and 
the block data USECDM in overlay (4, 0). 

Overlay (1, 0) - Entry STCN 

Overlay (1, 0) is the primary overlay for the program input and initial 
construction of the data tables and the unrefined calculation net. Program 
STCN serves as a driver routine to call in subordinate secondary overlays. 
These links in turn read the input and build the initial tables and calcu- 
lation net. Subroutines used concurrently by secondary links also reside 
in overlay (1, 0). 
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Overlay (1, 1) - Entry STCNR 


Overlay (1, 1) consists of the subroutines to read the card and/or 
tape input and perform preliminary processing. General input data are read 
in REDINP and boundary and channel data are read in RBD and RCD respectively. 

The channel data input table is also constructed in RCD. In the case of the 
boundary input, if angles are not specified with the coordinates, the SM9TH 
routines are called to calculate angles by fitting a beam through the specified 
coordinates. Also, the NACA Series 1 Cowl coordinates are stored internally 
and may be used to generate an "analytic" contour for a specified boundary 
segment. The input quantities to select these options are discussed in Section 
7.1, and on the input sheets (Section 12.2). 

Overlay (2, 1) - Entry STCN2 

The principal function of overlay (2, 1) is to build the boundary table 

and the leading edge/trailing edge table using the coordinate and angle data 

from the preceding link. This procedure is accomplished in subroutine BLDTAB. 

If specified by input, this overlay is recalled to initialize the matrix for 
the solution on the far field boundary (FRFDNZ) . 

Overlay (1, 3) - Entry STCN3 

The bulk of initial table construction is performed in subroutine BLDTBS. 
Specifically, the orthogonal channel table, the streamline table, the station 
table, the flow adjustment table, and the wake displacement thickness table are 
built in this routine. Subroutine BC0NV is called to build the convected 
properties table as an intermediate task. In the process of table construction, 
the values of Z, R, S 2 » and V for the unrefined grid are placed in the field 
tables . 

Overlay (2, 0) - Entry STCB 

Overlay (2, 0) is the primary link for calculation of the coefficient B and 
the RHS of the streamline correction equation using the "flow balance" relations 
The subroutine FL0BAL and TTPT are included in this link, since they are used 
by the output routines to calculate final velocities in the field prior to 
printing. 
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Overlay (2, l) - Entry STCS 

Overlay (2, l) contains the subroutines STALOO and BRHS to loop through 
the stations and calculate the coefficients B and RHS for the matrix solution 
of the streamline adjustment equation. Included also is the subroutine ADJWF 
which adjusts nozzle flow rates prior to entering the flow balance calculation. 

Overlay (2, 2) - Entry STCtfl 
Overlay (2, 3) - Entry STCW2 

These two overlays are loaded at the end of the problem solution to produce 
the printed output and generate a restart tape if requested. Subroutine WRIA 
prints the parameters defined by input, program storage utilization information, 
solution convergence level as well as generating the restart tape. Subroutines 
VfRieUT and WRIBDY print out the solution at the field points and the boundaries 
respectively. < 

Overlay Q, 0) - Entry STCXX 

Overlay (3, o) is loaded prior to the execution of the flow balance. The 
prime functions of this link are to adjust the streamline positions (ADJSL) , 
refine the grid (REFINE), calculate streamline curvatures (SLC), orthogonalize 
the grid (PTMGVE) , and calculate velocities on the far field boundary (FARFLD) . 
The far field routines are only called if these options are in effect. REFINE 
will not be called if grid refinement limits have been reached, the flow bal- 
ance error is too large, or the maximum field point limit has been met. 

Overlay (4, 0) — Entry STCM 

Overlay (4, 0) contains the subroutines MCOEF and IAD to set up and solve 
the matrix equations for the streamline adjustment. 
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if grid refinement 

and final tolerance additional grid refinement 

are satisfied if rough tolerance is satisfied 


,t> v 


THE STREAMTUBE CURVATURE PROGRAM 


Routines 


Function 



EXIT 


Read Input. 

Build Tables and 
determine first crude 
calculation net. 

Calculate distances between 
points along streamlines. 

Refine grid by adding new 

orthogonals and streamlines as required. 


Calculate streamline curvatures, 
angles, and distances between 
points. Locate stagnation points. 

Move points along streamlines 
to obtain orthogonal lines. 

Calculate velocities on the 
far field boundary. 

Adjust nozzle flow rates. 

Station loop to calculate 
coefficient B and RHS using 
"flow balance" equations. 

Calculate matrix coefficients 
Including far field B. C. 

Solve matrix equations 
for streamline movement. 

Adjust streamlines. 


Print output. 


Fig. 5-1 Program Flow Chart 


Output 
Z, R, S2, V 

SI 

Z, R of refined grid 

CURV, PHI1, SI 

Z, R, CURV, PHI1, SI, S 
UDN 

B, RHS, ES2 

Al, A2....A8 

DS2 

Z, R 
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Figure 5-2 STC Overlay Structure 




Fig. 5-3 Curve Fit Algorithm 
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Point Movement for Obtaining Lines 



Small Perturbation 
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Fig. 5-6 Illustration of Far Field and Near Field Solution Domain 



? 

6.0 PROGRAM NOMENCLATURE 

Communication between the subroutines in each overlay is accomplished by the 
use of labeled consnon. with few exceptions, the majority of labeled common 
storage areas are located in the main overlay. The size and data in these 
blocks are initialized by the block data subprograms USECDG and STCBLK. 

The principal labeled common blocks are given alphabetically in the following 
section. Labeled commons used by general purpose utility subroutines are not 
included in the tabulation. In most cases, the use of these commons is given 
in the/listing of the pertinent subroutine ( Volume III). 

Within each block, variables are listed according to the position occupied in 
the block. In some cases, the variable name may differ between routines and 
a typical name is given. Also, several areas are used primarily as erasable 
temporary storage and are denoted as such in the description of the labeled 
common. Inspection of the source listing should indicate the particular use 
of these areas by a given subroutine. The pertinent dimension and type infor- 
mation are included with the variable name (R = Real, I = Integer, L = Logical) 
Variables normally containing BCD data are typed as H = Hollerith, even though 
they may have real or integer names. 

v ■ ' 

r"' 
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BDU H 1 Upper boundary name associated with point 

NUSED - 491 Additional tables having the preceding format 
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7.0 STC PROGRAM INPUT/OUTPUT 

The following sections are concerned with the input to the STC program, 
special user instructions, and the output produced by the STC program. The 
standard system files INPUT and OUTPUT (TAPE6 = OUTPUT) are used for card 
input and printed output respectively. Additional data files designated 
as TAPE1 (input) and TAPE2 (output) may be used by the programs. In general, 
these files will reside on magnetic tape. 

7 . 1 Program Input 

Input sheets for the STC program are given in Section 12.0, along with 
special notes pertinent to the use of these sheets. Optional program input, 
not normally required, will be described in this section. Input data 
may be in the form of punched cards or a magnetic tape file (output file 
from a previous execution of the STC program) . Data read from a magnetic 
tape file may be selectively over-ridden or augmented by input cards. Four (4) 
distinct card input sets are read by the program and are: 


1.) 

Input 

sheet 

0 

Identification information 

2.) 

Input sheet 

1 

Overall input data 

3.) 

Input 

sheet 

2 

Boundary coordinates 

4.) 

Input 

sheet 

3 

Channel flow properties 


The first input set, consisting of identification information, is read once 
in a given run using fixed field format (6A10). The remaining input sets 
consist of a header card followed by a NAMELIST input list $A. Standard 
FORTRAN IV NAMELIST (Volume III) is used to read these latter lists. 
Successive cases may be run using only input sets 2.) through 4.). In all 
cases, the channel flow properties (input set 4) may or may not be present. 

The input parameters for each set are given in the following section. 
Included in these descriptions are the input items appearing on the input 
sheets as well as controls for special program options and the modification 
of preset tolerances. 



7.1.1 Identification Information Input Sheet 0 

The first three (3) cards of the input deck consist of the name and address 
of the user and the problem identification. 


Card No. 

Cols . 



1 

2-61 

User name (1-60 

alphanumeric characters) 

2 

2-61 

User address or 
characters) 

location (1-60 alphanumeric 

3 

2-61 

Problem identification (1-60 alphanumeric 
characters) 


Blank cards may be substituted for input quantities not required. 


7.1.2 Overall Input Data Input Sheet 1 

The first card Of this input set is a header card consisting of a 1 in 
column 2, the word STC starting in column 4, and a T or an F in both columns 
14 and 24. 


Card Column 


Description 


2 

4-6 

14 

24 


1 - Denotes overall data input 

STC - Denotes program name 
Input tape? (T or F) 

Output tape? (T or F) 


The header card is followed by the NAMELIST $A and the associated overall 
input data. The NAMELIST input is terminated with a $ in column 2. 


Input Parameters for General Usage 


Parameter 


Description 


Preset 

Value 


MACHO 

Mach Number 


TSO 

Ambient Temperature 

1.0 

PSO 

Ambient Pressure 

1.0 

RG 

Gas Constant 

1.0 

GAM 

Ratio of Specific Heats 

1.4 

RHL 

Highlight Radius 

1.0 
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Parameter 


Description 


Preset 

Value 


RM 

Maximum Body Radius 

1.0 

TTE 

Body Closure Tolerance 

0 . 

AXI 

Problem Type T = Axisymmetric 
F = Plannar 

T 

GR(1) 

Table of Transverse Coordinates for 
Grid Refinement Tables (up to 10 values) 


SGR(l) 

Table of Transverse Grid Refinement 
Criteria (up to 10 values) 

1 . 

NGR 

Number of Entries in the GR, SGR tables 

1 

GZ(1) 

Table of Axial Coordinate for grid 
Refinement Tables (up to 10 values) 


SGZ(l) 

Table of Axial Grid Refinement 
Criteria (up to 10 values) 


NGZ 

Number of entries in the GZ, SGZ tables. 

0 

VMG1 

Maximum Mach Number Increment Between 
Grid Points in Streamwise Direction 

0.1 

VMG2 

Maximum Mach Number Increment Between 
Grid Points in Normal Direction 

0.1 

MAXIT 

Maximum Number of Grid Refinements 



Optional Input 

The following input quantities are considered optional input in that 
they are not normally required for execution of the STC program. These 
items consist, in general, of input controls for special program options 
and input data to modify preset or initialized constants and parameters. 


Stagnation Properties, Optional 


Parameter 


Description 


Preset 

Value 


PTO Stagnation pressure 1.0 

TTO Stagnation temperature 1.0 

The total pressure and total temperature may be input by specifying PTO and 
TTO normalized by the free stream static temperature and pressure. Alternately 
the total conditions are calculated from the Mach number, MACHO, and the 
ambient pressure and temperature, PSO and TSO. 
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Special Controls for Supersonic Flow, Optional 


Parameter 


Description 


Preset 

Value 


SSFML Formula number for calculation of 1 

supersonic curvatures 

= 1 for 3 point parabola 

= 2 for 4 point piecewise cubic (i.e. beam) 

= -1 for 3 point parabolas for both 

subsonic and supersonic points 

SSEF Supersonic entering flow along the F 

upstream boundaries? T or F 

SSEANG Entering flow angle for SSEF * T (degrees) 0. 

SSDF Supersonic flow downstream of a choked F 

station? T or F 

SSFEND Supersonic beam downstream end .75 

condition (SSFML ■ 2) 

0. = Parabola 

1. ■ Cubic 

SSFND1 Supersonic beam upstream end .75 

condition (SSFML * 2) 

0. ** Parabola 

1. *> Cubic 

BRLX B - Relaxation factor 1. 

TSIC Number of points to be used for 2. 

transonic interpolation of curvature 


Boundary Conditions for Streamline End Points, Optional 
Parameter Description 


Preset 

Value 


NBC 


AFC 


NBC(l) = Upstream, NBC(2) = Downstream; 
NBC = 0 Angle specified 

NBC = 1 Curvature specified 

NBC = 2 Ratio of the rate of 

curvature at the end point 
to the value at the next 
end point is specified 


AFC(l) = Upstream, AFC(2) e Downstream; 0. 

Angle, curvature or ratio of rate 
of change in curvatures 
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Matrix Solution Controls, Optional 


Parameter 

IADM 


RHOBAS 

RHOAMP 


MAXSWP 

TOLRL 


Description 


Sweep parameter 

0 = IAD 

1 = Orthogonal block relaxation 
-1 = Streamline block relaxation 


Base acceleration factor 
parameter, p 

D 


Half amplitude of the acceleration 
factor, p. 

Note: A 


P 


P B + 


2p. sin 
A 


nn 

2*^NM 


Maximum number of sweeps 

Relative tolerance for matrix solution; 


|ADS2| < T0LRL*DS2 
1 1 — max . 


Flow Balance Solution Controls, Optional 


Parameter 


Description 


T0LES2 

TOLINR 

NINNER(l) 

CNVF(l) 

NODENS 


Relative tolerance on maximum point 
movement predicted by the flow balance 

Inner iteration tolerance on streamline 
movement 

Number of inner iterations (without 
grid refinement). Specify up to MAXIT 
values 

Fractional percentage of total point 
movement to be used at a given refinement 
level. Specify up to MAXIT values. 

Grid refinement level through which 
incompressible flow properties (y=“) 
will be assumed. 


Preset 

Value 

0 

0.5 

0.5 

200 

.001 


Preset 

Value 


16*10 

10 * 1.0 

0 
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Grid Refinement Parameters, Optional 

The following input items may be used to control the length of ortho- 
gonal lines and streamlines as the grid is refined. The use of these 
parameters is discussed in the following section on "special user instructions" 


Parameter 

Description 

Preset 

Value 

CRXSL 

Used for extending new streamlines 

.375 

CRXOL 

Used for extending new orthogonal 
lines across a subsonic region 

.375 

CRXSS 

Used for extending new orthogonal 
lines across a supersonic or mixed 
region 

.125 

CRXE 

Used for extending new orthogonal 
lines which cross an expanding 
flow sonic line 

0 

CRXC 

Used for extending new orthogonal 
lines which cross a supersonic to 
subsonic compression line 

0 

Optional Print 

and Control Options 


Optional diagnostic print and special logic controls are 

located in 

common /CPRINT/. Normally these variables are set to 0.. 

As indicated 


the following section, the print controls will provide output from a number 
of routines. The setting, pertinent subroutine, and resulting printout are 
indicated. Several variables control logic flow only. 

Parameter Description - 

PRTES2 Subroutine BRHS ’ 

=2 Print B,RHS,DS2,Z,R,PHI,CURV, and ES2 
for the two stations with maximum ES2. (DS2 
is the value computed during the previous 
iteration.) 

>2 Print the above information for 

PDUM(8) <C 1 lPRTES2. (B,RHS and DS2 are 
the values computed during the previous 
iteration.) 
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PRTA 

Subroutine IAD 

Print matrix solution for DS2 at 
end of each sweep. This indicator 
is turned on internally if 
NSWP > MAXSWP - vSm - 2 

PREFIN 

Subroutine REFINE 

Print diagnostic during streamline 
refinement for 
PDUM(8) <S 2 1 PREFIN 

PREFN2 

Subroutine REFINE 

^0 Print diagnostic during orthogonal 
line refinement for 
PDUM(8) <$ <_ PREFN2 

SSONIC 

Subroutine BRHS 

=2 Suppress area scaling for ES2 
calculation 

PDUM(l) 

Subroutine SLC 

>1 Print ZB , RB , ANGD , CURVD during 

iteration for stagnation point location 

=1 Print results of curvature calculation 

=2 Print results of supersonic curvature 
calculation 

=4 Print curvature results if point is a 
stagnation point or the termination of 
a partial orthogonal 

PDUM(2) 

Subroutine IAD 

Print final matrix solution for DS2 

PDUH(3) 

Subroutine MCOEF 

^0 Print influence coefficients A 2 »...Ag 

PDUM(5) 

Subroutine MCOEF 

>1 Print influence coefficients, G, relating 
streamline point movement to the negative 
of curvature 

=3 Print above information for partial 
orthogonal points only 

=4,5 Print above information for 1st and 2nd 
points upstream of leading edge 

PDUM(6) 

Subroutine ADJWF 


?^0 Print results of flow adjustment iteration 
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Description 


Parameter 

PDUM(7) 

PDUM(8-9) 

PDUM(IO) 

PDUM(ll) 

PDUM(12) 

PDUM(13 

PDUM(18) 

PDUM(20) 

PRPRN 


Subroutine REFINE 

1*0 Print velocities, stream-function after 
refinement 

Subroutine FLOBAL 

Print flow balance related data if 

PDUM(8) <£. <_ PDUM(9) and 

PDUM(8)*PDuM(9) + 0 

Subroutine BUI LPT 

1*0 for printout of tables described in 

Section 5.1 after they are first constructed 

-2 for printout of all tables at end of 
solution 

Subroutine FARFLD 

1*0 Print coordinates, slope, velocity on 
the far field boundary 

Subroutine MCOFF 

^0 Use spline influence coefficients even 
if SSFML - -1 

Subroutine FARFLD 

^0 Use "linear" interpolation on the far field 
boundary 

7*0 Use "least squares parabolic" interpolation 
on the far field boundary 

Subroutine FLOBAL 

=1 for calculation of normal shock total 

pressure loss for transonic compressions 

Subroutine FRFDNZ 

1*0 Print matrices used in far field solution 

Subroutine WRIOUT 

^-1 Print field information at each station 


7.1.3 Boundary Coordinates Input Sheet 2 


The first card of this input set is a header card containing a 2 in 
column 2, the word BDY starting in column 4, the boundary name in column 14, 
and the channel name in column 24. 
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Card Column 

Description 

2 

2 - Denotes boundary input 

4 

BDY - Denotes input type 

14 

Boundary name, 1-6 alphanumeric 
characters 

24 

Channel name, 1-6 alphanumeric 
characters 


The header card is followed by the NAMELIST $A and the associated input for 
the specified boundary. The NAMELIST is terminated with a $ in column 2. 


Parameter 


Description 


UPPER 


ZRONLY 


Boundary position 
UPPER - T 
UPPER - F 
Geometry indicator 
ZRONLY = T 
ZRONLY - F 


Upper boundary 
Lower boundary 

No surface angle input 
Surface angle input 


Two options are available for the input of the boundary geometry. Either 
the coordinates and the surface angle (measured from the positive x-axia) 
may be input, or the coordinates alone may be input. In either case the 
points must be input accurately. The first option (ZRONLY =■ F) is preferred 
With the second option (ZRONLY ■ T) , a beam is fit to the input points to 
determine the angles and, in the process, the beam fit angles and curvatures 
are printed. To determine whether a suitably smooth curve has been fitted 
to the points, the user should examine the curvatures and make sure they 
are reasonable. In general, the points should not be closely spaced except 
as required in regions of high curvature. In these regions, the angle 
change between points should be less than 25 degrees . 


The NACA Series 1 cowl coordinates are stored internally, and may be 
selected by specifying ZRONLY ■ T and B(l) = 991, 1, followed by the high- 
light coordinates and the coordinates of the maximum diameter. If the cowl 

i 

is to be extended beyond the end of the Series 1 contour, the "maximum 
diameter" coordinate is repeated and then other coordinates are listed. 
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Boundary coordinates are normally input in tabular form using the B block 
as specified on input sheet 2. All points must be listed in the streamwise 
direction for each boundary. Points at sharp corners must be listed twice, 
once for each angle which exists at the point. Normally, pressure and Mach 
number distributions will be printed at each orthogonal intersection with 
the boundary. Orthogonals may be forced to coincide with boundary points 
by listing the point twice in the input and setting DBLPTS » 0. (see optional 
boundary input). 

The column names in the B arrary are Z (or X), R (or Y) , and ANGD. 


Parameter Description 

B(l) Input block for boundary data 

column 1 Z or X coordinate 
column 2 R or Y coordinate 
column 3 Angle or slope of surface 

If desired, the data in the B array may be input in "free form" using the 
symbolic names associated with the appropriate columns of the B block. 

Parameter Description ' 


R 

Y 

Z 

X 

ANGD 


Normal (radial) coordinate (axisymmetrlc) 
Normal (vertical) coordinate (Planar) 
Axial coordinate (axisymmetrlc) 

Axial coordinate (planar) 

Angle or slope of surface 


Boundary data input via the B block will override corresponding data input 
in the "free form". 


Optional Input 

Optional input parameters, not shown on input sheet 2, are available to 
apply linear transformations to the input coordinates. Also, input control 
may be specified to force orthogonals to be placed at selected boundary input 
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points. This is accomplished by repeating the coordinates in the boundary 
input and specifying DBLPTS ■ 0. Normally, DBLPTS *» .01 forces removal of 
extra orthogonal stations where the angle discontinuity is less than .01 
degrees . 


Parameter 


Description 



DBLPTS 

ROTATE 

ZPIVOT 

RPIVOT 


Double point tolerance. Double points 0.01 

will be deleted if the angle discontinuity 
is < DBLPTS (degrees) 

Angular rotation in degrees 

Coordinates of the pivot point of 
rotation 


SCALE 

t 

ZTRANS 

RTRANS 


Multiplicative constant to be applied to 
the coordinate data 

Translation Increment in axial direction 

Translation increment in radial/vertical 
direction 


7.1. A Channel Flow Properties 


Input Sheet 3 


The first card of this input set is a header card containing a 3 in 
column 2, the word CHN starting in column 4, and the channel name starting 
in column 14. 


Card Column 


Description 


3 

4 
14 


3 - Denotes channel input 

CHN - Denotes input type 

Channel name, 1-6 alphanumeric 
characters 


Channel flow properties are specified in the following NAMELIST A, 
which is terminated with a $. 



Parameter 


Preset 

Value 


Description 


GAM 

RG 

TTO 

PTO 

MACHO 

TSO 

PSO 


Ratio of specific heats 
Gas constant 

Total temperature normalized by 
freee stream static temperature 

Total pressure normalized by free 
stream static pressure 

Mach number 

Static temperature normalized by 
free stream static temperature 

Static pressure normalized by free 
stream static pressure 


1.4 


1.0 


1.0 


1.0 


Specify either TTO and PTO (if known) or MACHO, TSO, and PSO if the Mach 
number is known. When these options are not used, the second option is 
employed with the free stream values of PTO and TTO computed from the 
MACHO supplied in the overall input data (sheet-1) . 


Parameter 


Description 


Preset 

Value 


AO Nondlmensional flow area normalized by 

the product of it and the square of the 

highlight dimension. For plane cases 

input AY /AY™ . For inlets AO I mass flow 
. o rli-> 
ratio. 


The flow is computed using the total properties as determined from the 
supplied Mach number, MACHO, and the flow area AO. If input data are not 
supplied for a given channel, the reference properties on input sheet-1 
will be employed with the area calculated at the entrance station. 


Individual sets of boundary and channel data need not be input in any 
particular order. For example, boundary input for the upper and lower 
surfaces bounding a given channel (input sheet 2) may be followed by the 
input for the pertinent channel (input sheet 3) . While arbitrary placement 
of these input sets may be used, it is well to develop some standardized 
conventions for the sequential input of boundaries and channels. This will 
be discussed further in the following section on special user instructions. 
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7.2 Program Output 

The output from the STC program may be logically divided into the 
following six sections: 

1. ) Card input and preliminary printout. 

2. ) Input and calculated boundary coordinates and angles. 

3. ) Solution history 

4. ) General input and output data 

5. ) Flow field data along orthogonal lines 

6. ) Calculated flow data along field boundaries and final channel 

momentum balances . 

The above sections appear sequentially in the output except when 
PRPRN = -1. In this case, the print of the flow field data along orthogonal 
lines is eliminated and the calculated boundary information follows the 
general input/output section. 

7.2.1 Card Input and Preliminary Printout 

\ 

The initial section of output consists of a card image print of the 

problem input and a designation of the tape Input/output file selections; 

T T 

viz, TAPIN = and TAPOT = _ . Upon completion of the card image print, 

f' F F 

the file TAPE5 is rewound to its original position. 

7.2.2 Input & Calculated Boundary Coordinates & Angles 

Boundary coordinates are processed and printed as they are read. The 
columnar output displayed for each boundary consists of the input coordinates 
<X,Z), (Y,R) and the input or calculated body slope (ANGD) in degrees. 
Preceding these items is a bold heading specifying BOUNDARY COORDINATES, the 
boundary name (BDY = name), the adjacent channel (CHN = name), and whether 
the surface is an upper (UPPER = T) or a lower (UPPER ■ F) boundary. When 
the ZRQNLY = T option is in effect, intermediate printout will be produced 
defining the results of the beam curve fit to the input coordinates. Since 
no smoothing is applied to these coordinates, the only meaningful output is 
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the (X,Z) (Y,R) and the calculated angles (ANGD) . Break points in boundaries 
(double points) enable the curve fit of a boundary in segments. The consoli- 
dated output, consisting of the collated boundary data for all segments, 
follows the above intermediate printout. The NACA Series 1 cowl coordinates 
(ZRONLY = T) appear as 40 points in a form similar to that described previously 

Normally, the solution history (described in Section 7.2.3) is printed 
as the next section of output. When far-field boundary conditions are applied, 
however, a comment is inserted in the output designating the original trans- 
verse location and streamwise extent of the far field boundary. Prior to 
calculation, the upstream and downstream limits of the far field streamline 
are extended to insure that the velocities are well behaved at the stream- 
wise ends of the STC integration region. The coordinates of the resulting 
far field boundary are printed below the above comment. 

7.2.3 Solution History 

The solution history output provides a convenient summary of the problem 
history as the solution proceeds through major grid refinements and inner 
iterative improvements of the flow balance. Included in this section are the 
parameters defining the grid refinement, the inner iterations , the grid 
othogonallzatlon, the flow balance errors, and the matrix solution. The 
pertinent variables in their literal order of appearance in the solution 
history list are: 

Description 

Grid refinement level 

Number of field points 

Inner (flow balance) iteration level 

Convergence factor on streamline movement 

Root mean square point movement along streamline 
before damping 

Maximum point movement along streamline before 
damping 

Maximum point movement along streamline after 
damping 


Variable 

NREFIN 
GRID PTS 
INRCTR 
CNVF 
RMS-DSl 

MAX-DS1 

MAX-DS1 
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Variable < 
% 

LIM-ES2 

MAX-ES2 
MAX-DS2 l 

N SWEEPS 


_ Description - 

Limiting streamline ES2 for satisfactory 
solution of the flow balance 

Maximum ES2 predicted by flow balance 

Maximum streamline correction predicted by 

the matrix solution for DS2. =0. when INRCTR = 0 

Number of iterative sweeps for solution of the 
matrix equations . = 0 when INRCTR => 0 


7.2.4 General Input and Output Data 

/, 

r i 

The output in this section consists of general input and a summary of ... 
selected output parameters. All of these items have been defined previously 
in Sections 7.1 and 6.0, hence only the general subdivisions of this 
printout will be listed. 


JX 



a. ) General Input 

b. ) Streamline End Conditions 

c. ) Supersonic Curvature 
, Parameters 

d. ) SubsoniJ/Supersonic 

Branch Selection 

e. ) Grid Size Criteria 

f. ) Memory Utilization 


Section 7.1 
Section 6.0 
Section 6.0 

Section 6.0 

Section 7.1 


The memory usage in terms of field points, STC table storage, and 
number of streamlines is compared with the maximum available values 
of these parameters. 


g. ) Convergence Data Section 7.1 

h. ) Special Boundary Options 

The GE version of STC contains a number of special boundary options 
in addition to the hard wall and the far field. Only the hard wall 
and far field have been explicitly implemented in the NASA Langley 
version of STC, due to limitations on the hollerith input capabilities 
of CDC NAMELIST input. 

i. ) Matrix Solution Section 6.0 
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j. ) Contents of Channel Section 6.0 

Input Table 

k. ) Channel Flow Rates, 

Pressures and Temperatures 


Included in this output are the channel flow rates and the pressures 
and temperatures in each individual channel. Note that the flow rate in a 
given channel may be adjusted to satisfy the Kutta condition at a trailing 
edge. When dimensional properties are input, the columns denoted PT/PSO and 
TT/TSO have units compatible with the dimensions of the input gas constant 
RG. 

7.2.5 Flow Field Data Along Orthogonal Lines 

Flow field data on orthogonal lines are printed when PRPRN ^ -1. The 
first portion of output on each page consists of the program name (STREAMTUBE 
CURVATURE PROGRAM) and the problem identification. The STC program uses an 
auxiliary coordinate system (£^, as a numbering scheme for orthogonal 
lines and streamlines. The coordinates = XII and “ XI 2 apply to the 
orthogonal lines and streamlines respectively. Following the problem identi- 
fication is a heading with the constant station value • XII along with 
the channels through which the given orthogonal passes. Primary (initial) 
orthogonals are flagged with ** after £. At the extreme right of this line 

A • 

of output is the flow type; viz., SUB ■ subsonic, SUP * supersonic, or 
CHOKE = choked. Pertinent data from the field arrays are printed in column 
format. Stagnation points do not appear in the flow field output. 


Variable 

Description 

XI 2 

Streamline numbering coordinate £ 2 

STEM FNCT 

Dimensionless stream function; fraction of flow 
in channel 

X, Z 

Strearawise coordinate 

Y, R 

Transverse or cross stream coordinate 

PHI 

Flow angle, deg. 
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Variable Description 


CURV 

Streamline curvature 



PS/PO 

Ratio of local static 
static pressure 

pressure to 

reference 

PS/PT 

Local static to total 

pressure ratio 

TS/TT 

Local static to total 

b < 

temperature 

ratio 

CP 

Pressure coefficient 




P - P 
s 00 

1 “ J2 

-zyl? M 

/ ' 00 00 


Compressible 


P - P 
s 00 


P M 

00 OO 


Constant density 

( Y = 0) | 


MACHA < .1, 


tr 


MACH 

AREA 

PTQPTO 


Mach number 
Flow area 

Channel total pressure ratio 


For primary orthogonals, the axial and normal momentum fluxes are 
printed following the field data. These items are: 

Axial Momentum Flux = J V cos <}> dw + J ” (P-P g ) cos <(1 dA 


Normal Momentum Flux 




V sin 4> dw + 


J * (P-P g ) sin <(> dA 
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7.2.6 Calculated Flow Data Along Field Boundaries -Final Channel Momentum 
Balances 

The final section of standard output includes the flow parameters on the 
upper and lower streamlines bounding each channel. Normally the boundary 
data follows the field data. In the case of PRPRN * -1, however, the boundary 
information is printed immediately after the general input/output section. 

Three lines of identification information head the boundary data, and 
consist of the program title (STREAMTUBE CURVATURE PROGRAM), the specific 
case identification, and a designation as to whether the boundary is an upper 
or lower portion of the given channel. The streamline coordinate (XI2) is 
given following the channel name. Boundary flow parameters are printed at 
orthogonal intersection and consist of the following items: 

Description 

Orthogonal (station) numbering coordinate (£^) 
Distance along boundary 
Axial coordinate 
Vertical or radial coordinate 
Flow angle or surface slope 
Streamline or surface curvature 
Local static to reference pressure ratio 
Pressure coefficient 

C - 2(P-P )/yP M 2 
p s s o 

C = Rg(P-P )/P M 2 (Y = 0, constant 

^ 8 8 o density) 

C = 0 (M < .1 

p o 

Local static to total pressure ratio 
Local Mach number 

Ratio of local total to reference static 
pressure 


Variable 

XII 
S1W 
XW, ZW 
YW, RW 
AN GW 
CURVW 
PS/PO 
CP 


PS/PT 

MACH 

PT/PO 
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The cumulative drag/ thrust is printed in the column labeled CDPI. The 
force is normalized by the free stream dynamic pressure (q) and the maximum 
area based on RM and is given by: 


CDPI 


(p - p ref ) 


dA 


max 


The adjacent column (A-AMAX) /AMAX is projected area normalized by the 
maximum area. Following the column printout is the ratio of boundary 
total temperature to ambient total temperature. When the boundary is 
the approach stagnation streamline, the inlet additive drag is printed 
below the total temperature ratio. This force is also normalized by q • 
and AMAX and is printed for both the upper and lower stagnation stream- 
line. 

Momentum Balance 

The STC program evaluates the thrust/drag on each boundary surface bounding 
a given channel and then verifies these forces by performing an "overall" 
momentum balance for each of the fluid streams. The integral momentum 
balance output follows the boundary output for a given channel. This 
information consists of the entering axial momentum, the integrated 
pressure forces on the boundaries, and the leaving axial momentum. The 
discrepancy between the leaving momentum and the sum of the inlet momentum 
and the pressure-area forces represents the net error in the calculation. 
This error may be attributed to inaccuracies in the computed pressure 
distributions or, perhaps, to insufficient refinement of the calculation 
grid for adequate resolution. 
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8.0 PROGRAM MESSAGES AND ERROR CODES 


8,1 Parameters Describing Solution Status 

The solution history output, described briefly in Section 7.2.3, contains a 
summary of how the solution is proceeding. By reviewing this solution status, 
the level of grid refinement and the degree of solution convergence can be 
readily determined. In brief, the solution history is like a "hospital 
patient's history" in that the general health of the calculation procedure 
is charted here. (See Section 11 - Sample Cases). 

i 

Ihe grid r refinement level is specified by the number of times the grid has 
been refined (NREFIN) and by the number of GRID PTS in the flow field. The 
number of grid refinements is controlled by MAXIT. The number of points in 
the flow field is dependent on the grid refinement criteria specified and 
MAXIT. 

The number of inner iterations at any grid refinement level is shown by 
INRCTR. These are controlled by built-in tolerances on the flow balance, 
but the internal program control can be overridden by inputting NINNER(m) 
where m refers to the level of grid refinement. 

The system of built-in tolerances is shown in Figure 8-1. The iteration 
tolerances are controlled by 70LINR and by T0LES2 which are input. The 
system starts with the maximum streamline movement (MAX-ES2) demanded by the 
"flow balance" in FLOBAL and compares this to a specified limit, L1M-ES2. 

If the MAX-ES2 is within limits, the "flow balance" is satisfied and more 
grid refinement is called for (if the limit has not been reached); otherwise, 
inner iterations continue until the specified tolerance, LIM-ES2, is met 
(or the specified NINNER limit is met). The specified tolerance, L1M-ES2, 
is controlled by TOLINR and the current grid size. Thus, the "flow balance" 
tolerance becomes smaller as the grid refinement proceeds. 

If the grid refinement (MAXIT) criteria is satisfied or the field point limit 
is reached, a second tolerance (T0LES2) defines a final series of inner 
iterations. This permits a looser tolerance (TOLINR) during the flow field 
grid development and then a final tight tolerance (T0LES2) for actual flow 
field definition. The summary of this iteration tolerance is shown as NREFIN, 
GRID PTS, LIM-ES2, and MAX-ES2. A scan of the two columns in the flow balance 
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(See Section 11 - Sample Cases) shows the decreasing values of MAX-ES2 as the 
solution proceeds. Any erratic behavior in MAX-ES2 indicates convergence 
problems (the patient's heart ia behaving poorly). 

During an inner iteration, the streamline displacement (or correction) 
equation is solved to adjust the streamlines in the flow field. The results 
of the matrix solution are shown as a maximum streamline adjustment (MAX-DS2) 
which is actually made and the number of sweeps (NSWEEPS) through the matrix 
solution (IAD or LRELAX) to reach the final solution. Also shown is the cumu 
lative time to complete the calculations at each level of refinement. 

Three other columns, titled "orthogonolization", summarize the orthogonal ad- 
justment in PTMOVE. This history is informative, but does not identify so- 
lution problems usually. 

Under certain output options, the diagnostic printout includes a line of so- 
lution history at the end. Since the history is not summarized neatly on 
one page, care must be exercised in seeking out the solution history inter- 
spersed with the diagnostic output. 
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STC ITERATION TOLERANCE 



2. Continue grid refinement if: , 

(a) MAX-ES2 < LIM-ES2 I 

or INRCTR > NINNER (NREFIN) 1 

and (b) NREFIN < MAXIT I 

and (c) grid refinement criteria 

has not been satisfied or 
field point limit has not 

been reached. • 

I 


KLOIiAL 


Otherwise I 
calculation | 
is complete | 


J 



LIMES2 = T$LINR* [current Grid Size) 

CLEN = characteristic grid interval 
length; it is equal to the a 
average of the SGR's & SGZ's. 


Continue inner iterations if: 

(a) MAX-ES2 > LIM-ES2 and 

INRCIR < NINNER (NREFIN) 

r (b) grid refinement criteria is 
satisfied or field point 
limit is reached and 
MAX-ES2 > T0LES2*CLEN 

T0LES2*CLEN = final tolerance 



MAXIT, NINNER, T$LINR & 
T0LES2 are input items. 


Fig. 8-1 Iteration Tolerance 
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8.2 STC Error Processing 


In general, the occurrence of an error during processing will cause 
termination of the problem. Each of the main links contain an error 
print routine which will print the STC tables and the field tables 
to assist in diagnosing the problem. When the error occurs in overlays 
subordinate to (0,0), an illegal computed GO TO statement is executed. 
This will cause a MODE 4 abort and a printout of the sequence of sub- 
routine calls preceding the error. 

8.2.1 Description of Table Print Format 


In general, the information stored in the STC tables consists of mixed 
BCD, floating point, integer, and logical data. A special routine 
(TABPRT) is used to print the arrays which may have a variable data 


format. The following conventions 
format on the CDC machines: 

Data Type 

Integer 

Floating Point 
BCD 

Logical 


have been established for the output 


Print Format 

Integer (112) Format 

E or F Format, depending 
on the magnitude of the 
numbe r . 

Hollerith characters 
A6 Format 

TRUE. Variables are 
printed as an integer -0 

FALSE. Variables are 
printed as an integer 0 


Note that a floating point 0, an integer 0, and a logical. FALSE, are 
all printed as an integer 0 since the data type, in this instance, is 
indistinguishable. The program listings ( Volume III) may be consulted 
to determine the true type of the data should this situation arise. 

A junk word BITS = 10^ is utilized in the STC program for testing 
purposes. This quantity is printed using OCTAL format (04) and appears 


9*2 


as 0000 . 



Since the data tables used in each overlay differs, the error print may 
differ slightly between links. The normal order of printed output consists 
of the general labeled commons, the STC tables (BDYTAB, STATAB, etc), the 
JMS table, and the field tables such as R, Z, VM, etc. The error print is 
best illustrated by an example of the typical output after an error in 
the flow balance (overlays (2,0), (2^,1)). The variables in the printout 
may be identified by referring to Section 6.0 or directly to the program 
listings. Referring to the example, the initial output consists of the 
data in labeled commons / ALLCOM/ and /CFB/. This latter common block is used 
as working storage by the flow balance routines. The STC tables BDYTAB, 

ADJWF and STATAB appear next. These data are packed in common/CHDATA/ , 
with the location designated by the limits set in common/ IXORIG/ . Note 
that the relative location within /CHDATA/ is printed as the extreme 
left column of numbers. Following the STC tables are the /SLTAB/ 
streamline table and the JMS table denoting the upstream /downstream print 
connections through the field. The JMS table is followed by the data stored 
in the field tables SI, S2, Z, R, PHI1, CURV, VM and B. The final printout 

•-ii.ij, 1 - ’ 

consists of the labeled common /ERASE2/, used as temporary storage by the 
flow balance routines. Following the output, a mode A error is initiated 
to obtain the subroutine backtrace. 

Common /TROUBL/ contains two logical variables ERR and INERR which may be 
set to .TRUE, during processing. In this situation, the job is terminated 
normally in the main program. The table printout is preceded by the 
comment ERR = T, ERRLOC = n^ , LENTRY = In this comment , n^ and 

are integers denoting the location of the last overlay call in the MAIN 
program and the current LENTRY setting (see listings). 

8.2.2 Location & Explanation of Specific Error Comments 

As mentioned previously, input or calculation errors will result in the 
termination of the problem. The locations (subroutine name) and description 
of error conditions are given in the following table. Included are suggested 
corrective measures (where possible), description of comments preceding the 
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table print, and a designation of whether the error is FATAL or NON- 
FATAL. One specific problem deserves further comment. When the 
iterative solution of the matrix equations for streamline point movement 
is in trouble (NSWP ■+ MAXSWP), a print of the solution matrix (DS2) and 
the maximum DS2 will start appearing in the printout. Subroutine IAD 
is structured to detect this problem and provide the iteration history 
over a complete cycle as the solution proceeds from over-relaxation to 
under- relaxation. In this case, the problem may usually be detected by 
inspection of the solution matrix. Appropriate corrections would be 
to input a RHOBAS . larger than the preset value or preferentially sweep 
the field in a single direction. 


Subroutine Location Type Specific - Description and/or Corrective Measures 

Name Nearest Comment , 

Statement No. 
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ISBOT 510 F No 

TTPT 125 F No Properties for requested channel not present in convected 

property table. 
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Subroutine Location Type Specific Description and/or Corrective Measures 

Name Nearest Comment 

Statement No. 
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Subroutine Location Ty pe Specific Description and/or Corrective Measures 

Name Nearest Comment 

Statement No. 
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9.0 OPERATING PROCEDURES 

The STC program described herein may be run on any Control Data 6400/ 

6600 machine operating under SCOPE 3.0 or a higher level operating 
system. In general, operating procedures and control card set-ups will 
differ from site to site. The following comments on program modifications, 
deck set-ups, and operating instructions are restricted to the program as 
installed at the NASA Langley Research Center. Minimal changes should be 
necessary for successful Installation at other CDC sites. 

9.1 General Operating Procedures 

Two versions of STC are available at the LRC computer site, both residing in 
source and absolute binary form on permanent file (data cell) disc storage. 

The first version has a maximum flow field data storage of 300 grid points, and 
may be run from the LRC remote input stations using a job field size of 
70gK locations. It is expected that this version will be used to run 
problems through the first 2 or 3 grid refinements to assure that the 
program is operating satisfactorily and the input is correct. For problems 
of nominal size, it will be necessary to run a larger version of STC at the 
central site (field length < 70gK) . 

The larger version of STC may be created by recompilation (FORTRAN IV RUN 
compiler) of the two block data subprograms which set the size of the STC 
table and field point storage (USECDG) and the size of the influence 
coefficient arrays used in the solution for the streamline adjustment 
(USECDM) . As a guide for setting the size of these labeled commons , a 
bench-mark field length of 100gK storage is adequate for a problem of 
730 grid points. The commons and arrays which must be altered to increase 
the problem size are: 
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1.) 

Block 

Data USECDG 

(overlay 0,0) 

Common 

Variable 

Table/Description 

/CDS2/ 

C12 


Field Table - Streamline adjustment 

/CRHS/ 

C16 


Field Table - Right hand side of matrix 
equation. 

/CHDATA/ 

C9 


STC Table Storage 

/CCURV/ 

Cll 


Field Table - Streamline Curvatures 

/CPHI1/ 

C14 


Field Table - Flow Angle 

/CS1/ 

C17 


Field Table - Distance along streamlines 

/CS2/ 

C18 


Field Table - Distance along orthogonal 
lines 

/SLTAB/ 

C8 


Streamline Table - Dimension 3* maximum 
number of streamlines 

/CM/ 

i 

1 

C13 


Field Table - JMS (mesh point connection) 
array 

/CB/ 

CIO 


Field Table - Coefficient B in matrix equation. 

/CZ/ 

C20 


Field Table - Axial coordinate (X,Z) 

/CR/ 

C15 


Field Table - Normal coordinate (Y,R) 

/CVM/ 

C19 


Field Table - Velocity 

As indicated in 

the preceding 

list, 11 variables comprise the flow field 


tables . These arrays should be dimensioned to the maximum number of mesh 
points. The exception is /CDS2/, which must be set to a minimum of 900 
locations. 

Note also that /CZ/ and /CR/ must reside next to each other in memory, as 
these quantities are used to determine the maximum limit for the number 
of field points. For problems with a large quantity of boundary input and/or 
a large number of anticipated orthogonal stations, it may be necessary to 
increase the size of the STC Table storage common /CHDATA/ beyond the 2046 
locations currently allocated. The streamline storage /SLTAB/ is set to 
accomodate 128 streamlines; this should be of adequate size for most problems. 
As in the case of /CZ/ and /CR/, the commons /CHDATA/ and /CEND/ must be 
located adjacent to one another in memory. These items are used to set 
the maximum limit for the STC table storage. 
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2.) 

Block Data USECDM 

(Overlay 4,0) 

Common 

Variable 

Description 

/CA2/ 

C42 

Curvature influence coefficient 
A2 

/CA3/ 

C43 

Curvature influence coefficient 
A3 

/CA4/ 

C44 

Streamline flow difference influ- 
ence coefficient A4 

/CA5/ 

C45 

Curvature influence coefficient 
A5 

/CA6/ 

C46 

Curvature influence coefficient 
A6 

/CA7 / 

C47 

Streamline flow difference influence 
coefficient A7 

/ CA8/ 

C48 

Streamline flow difference influence 
coefficient A8 


The dimension of the arrays in these commons should be the same as the 
dimension of the flow field arrays. 

9.2 Recommended Procedures for Program Maintenance & Modification 

As mentioned previously, a source and absolute binary copy of the of the 
STC program are stored on the permanent disc file (data cell) at the LRC com 
puter center. It is recommended that relocatable binary copies of each 
version be maintained on tape. As changes are made to the program, only 
those routines which differ need be recompiled. The relocatable binary 
tape may be updated using standard COPYL and COPYN functions. Subsequent 
execution from this file will result in the creation of a new absolute 
binary file which may replace the one currently on the data cell. The 
updated source decks may also be updated at this time. 

The source copy of STC contains *DECK cards as the first card of each 
subroutine. In running at other CDC sites, the source file may be 
used directly to initialize an update tape. 
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9.3 Standard Execution Deck Setups 

Normally, the absolute binary (overlay) program will be reloaded from the 
data cell for execution. The suggested deck set-up is: 

1 

(Job Card) 

(user Card) 

FETCH (A3727, address on data cell, BINARY, STC) 

COPYBR (INPUT, TAPES) 

STC. 

6 7g (EOR) 

[Input Data] 

6 7f , 9 (EOF) 

As indicated previously, this deck set-up is applicable at the NASA 
Langley Research Center computer site. At other installations, the 
absolute binary file may be reloaded from tape. For example, a typical 
deck set-up for use on the CDC Cybernet System is: 

1 

(Job Card) 

LABEL, STCOVLY, R, VSN = UT611, L ■= CE1. 

COPYBF (STCOVLY, STC) 

REWIND (STC) 

COPYBF (INPUT, TAPES) 

STC. 

6 7g (EOR) 

[Input Data] 

67 ^ (EOF) 

9.4 Restart of STC Program from Data Tape 

A partially completed STC problem may be restarted using the output data 
tape from the previous STC execution. The restart tape is obtained by 
specifying TAPIN as T on the general input header card and using a file 
card such as REQUEST or LABEL to assign IAPE2 to a given physical device. 
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Of course, the logical file TAPE2 must be assigned with "write permit" RING 
option. In the subsequent restart task, TAPIN must be specified as T on the 
general input header card. Also, the input file TAPE1 must be defined using 
an appropriate file card (NORING). The input for the restart case should 
consist of the identification cards, the general input header card, and the 
first $A namelist. For example, to restart a job after 3 major grid 
refinements and run to 5 grid refinements, the example input might 
appear as: 

2 

NAME =J. Smith 
ADDRES = LRC 

t 

IDENT = NASA 8 INLET 
2 A 14 24 

1 STC T F 

$A 

MAXIT = 5, 

$ 

Intermediate output at a given MAXIT level may be obtained by running a 
restart case and not assigning TAPE2 and TAPE1. In this instance, these 
files are assigned to disc and switched at the end of each total set of 
input. This option is also useful for changing program tolerances at a 
given stage of the calculation. 

The preceding example could be run as one case using the following input: 

2 

NAME = J. Smith 

ADDRES = LRC 

IDENT =NASA 8 1 INLET 

24 14 24 

1 STC FT 

$A 

MAXIT =3 

$ 

remainder of case input 
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1 STC T F 

$A 

MAXIT « 5, 

$ 

Successive restart cases may be stacked in sequence. On all but the 
final one, use 1 STC T T as the header card. 
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10.0 HELPFUL USER HINTS 


10.1 General Comments 

The following section is a compilation of information for the user. Many 
items have been covered elsewhere, but they will be repeated in this section 
for emphasis. The intent is to identify key input parameters so that com- 
puter solutions can be accomplished with minimum user problems or errors. 

The primary items of importance are: 

1. Smooth input geometry 

2. Grid refinement definition 
3- Iteration tolerances 

4. How to get out of trouble when the solution doesn't converge. 

10.2 Smooth Input Geometry 

The Streamtube Curvature method uses the boundary or body surface curvature 
as the geometry parameter which causes velocity and/or pressure gradients. 
Thus, the computer program is very sensitive to changes in curvature. Hence, 
an accurate definition of the boundary coordinates and local angles is imper- 
ative if surface curvatures are to be smooth and continuous. 

There are two methods to insure smooth input geometry; a) use an analytic 
boundary or function with continuous first (angle) and second derivatives 
(curvature) or b) smooth specified boundary coordinates by fitting a curve 
to the surface points. The definition of an analytic boundary is not always 
possible and, in some cases, the intersection of two analytic functions will 
have discontinuous second derivatives. (An ASME nozzle is an ellipse followed 
by a straight line. The discontinuous second derivative causes a pressure 
"blip"). Hence, a curve fit technique may be necessary in most cases. 

The best curve fit to use is a piecewise smooth cubic since the curvatures 
are evaluated in STC using the same type of curve fit. Hence, the definition 
of the coordinates and angles should come from some such geometry smoothing 
method. As part of the geometry definition, the curvatures should be checked 
for discontinuities or erratic behavior. 
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10.3 Grid Refinement 

The grid refinement criteria is explained in detail in Section 12. Some 
comments on use are appropriate with reference to Figure 12-3. 

First of all, the purpose of the grid refinement criteria is to maximize the 
grid refinement in areas of interest and to minimize the flow field grid in 
areas of smooth flow. By forcing flow field grid points to be used efficiently, 
computer table size and storage space is minimized and computer time is reduced. 
Also, the flow field grid can develop without over-constraining the streamline. 

The best technique is to first identify the region or regions in the flow 
field where a high degree of resolution is necessary. The geometric grid re- 
finement criteria in this area is set at a value equal to the grid size de- 
sired. As shown in Figure 12-3, a rectangular flow field section on the inlet 
lip is defined with SGZ = SGR. If a subsidiary plot of SGR vs GR and SGZ. vs 
GZ are set up, the input table is readily developed. In the radial direction, 
SGR is specified at the centerline, decreases linearly to the region of interest, 
increases linearly to the outer region value and is constant from there out. 

In the axial direction, a similar plot can be drawn. The coordinates of these 
two plots become the input tables for GZ/SGZ and GR/SGR. 

In some problems, the geometric grid refinement must be redefined if solution 
difficulties occur. A demand for excessive grid refinement in a local area 
can over-constrain the development of the flow field. Only experience will 
allow the user to make the optimum choice of grid refinement criteria for 
different problems. 

10.4 Iteration Tolerances 

lhe iteration tolerances, T0LINR and T0LES2, have been preset at recommended 
values for most problems. As explained in Section 8.1, these tolerances 
control the iteration logic for the grid refinement and the "flow balance" 
or inner loop solution. 

Usually, if the input geometry is properly defined and no errors are made in 
defining the flow properties (such as setting VARY = T in all the channels), 
the solution proceeds through several grid refinements with no problems. 
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When the grid refinement criteria and the iteration tolerances are incompat- 
ible, i.e., the tolerances are very small and the grid refinement has high 
resolution in an area of large streamline curvature, the solution will show 
instability in the definition of streamline adjustment (MAX-ES2 will not de- 
crease continuously). The remedy is to change the tolerances or the grid 
refinement or both. 

One good procedure is to keep the iteration tolerances relatively large during 
the flow field grid development and then restart the last grid refinement 
with smaller tolerances. Again, experience in use of the computer analysis 
will be necessary to provide guidance in how to do this. 

10.5 Trouble Shooting 

The Streamtube Curvature Analysis has been set up as a user oriented com- 
puter program for the solution of relatively complex transonic flow field 
problems. Many of the control variables have been preset to the most generally 
useful value. Hence, a majority of problems will be solved quickly and ac- 
curately if the inputs (geometry and channel flow properties) have been de- 
fined correctly. But, there will be instances when numerical difficulties 
will occur. 

Invariably, a large percentage of the numerical difficulties can be related 
to errors in defining the geometry or bad judgement in selecting the grid 
refinement or iteration tolerances. However, when large supersonic bubbles 
are present in a transonic flow, convergence problems in the inner loop have 
been experienced. A dump, very similar to that shown in Section 8 . 3 , will re- 
sult. The solution history gives the first diagnosis in that the level of 
grid refinement, the number of inner iterations during the last refinement, 
and the behavior of MAX-ES2 can be quickly reviewed. Often MAX-EST will 
start to converge and then show instability. By identifying the appropriate 
NINNER, the solution can be forced to the next level of grid refinement by 
rerunning the problem with the NINNER control inserted. 


123 


One good procedure is to rerun the program with the MOtIT set to the previous 
level of successful grid refinement and convergence, to get a complete print- 
out of the flow field, and then to restart the computer program with new iter- 
ation tolerances and/or grid refinement controls. The example cases in 
Section 11 show a problem solution which has been restarted with changes in 
iteration tolerances. 
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11.0 SAMPLE CASE - PROGRAM OUTPUT 

The following sample case is the analysis of an inlet lip at a free stream 
Mach number of 0.80. The geometry is defined as a straight centerline for 
the lower boundary of the inlet flow (BDY CNILN V2 and UPPER = F), the 
coordinates and angles for the inner surface of the inlet lip, the call for 
a NACA I external cowl lip (CLEX), and a farfield (FF) upper boundary. The 
grid refinement has been redefined as shown to eliminate some numerical 
problems. Hie iteration tolerance is T0LES2 =1.0. The output print option 
is initially set for a summary printout of boundary information only (PRPRN = -l). 
The grid refinement level is set as MAXIT = 8. The input is printed out as 
the first pages of the output. 

After eight grid refinements, the solution history and boundary flow con- 
ditions are printed out. The solution is restarted with MAXIT t= 12 and again 
the same printout is requested. 

The final iteration keeps the same grid refinement, but decreases the tolerance 
T0LES2 = .OOl for flow balance convergence. The print option is changed to 
get the full flow field. These results are printed out as shown. 

A comparison with test data is shown in Figure 11-1 for the solution after 
eight grid refinements (approximately 700 grid points) and after twelve grid 
refinements with the converged flow balance (1100 grid points). Ihe additional 
resolution near the leading edge of the external lip surface is evident. This 
is an inviscid solution so that the oscillations in pressure are defined. A 
boundary layer solution indicates local separation and the test data agrees 
with this. Once the inviscid solution gets downstream of the separation 
region, it matches the test data very well. 
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SUM-VM*COS(PHl)*orLOW « 73.70 SUM-VM*SIN (PHI > *DFLO* = 4.31 
SUM_(p_pso>«COS(RHn»OA = 76.92 SUM- (p-pSO) *SIN (pHl ) «DA = 4.50 
TOT AXIAL MOMENTUM FLUX = 150.62 TOTAL Y-MOMENTU* FLUX = 8.81 
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STREAMTUBE CURVATURE program 
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5UM-VM*C0S (PHI ) *OFLOW S 10025.1a SUM-VM*SIN(PHI > «DFLOW = 43.85 
SU"- (P-PSO) *COS <PHI ) *DA = -26.14 SUM- (P-PSOl *SIN (PHI ) *0A = -.11 
TOT AXIAL MOMENTUM ElUX s 9999.00 TOTAL Y-MOMENTUM FLUX * 43.74 
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-Q PROGRAM INPUT SHRPTR 



STREAMTUBE CURVATURE PROGRAM 

Aug. 1972 Page of 

Overall Input Data STC/Sheet-1 



input tape? output tape? 
T or F T or F 


1 STC 


Mach number, ambient pressure and temperature, fluid properties 


$A (1.) (1.) (1.) (1.4) 

MACH0= , TS0- , PS0= , RG= , GAM= 

Highlight radius, maximum body radius, body closure tolerance 

RHL= , RM= , TTE= , 

axlsymmetric or planar? 

(T) or F 

AXI= , 

spacial grid refinement criteria, see notes 

GR(1)= , » _» » * > _ 

SGR(1)= , , t » » » _ 

NGR= , 

GZ ( 1 ) = , ___ , 9 » > » _ 

SGZ ( 1 ) = , , - » » * » _ 


NGZ= , 

maximum Mach number increment between grid points 

normal 
direction 
( 0 . 1 ) 

VMG1= , VMG2= , 

\ I 

maximum number of refinements 
MAXIT= . ' : 


streamwise 

direction 

( 0 . 1 ) 




> 
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Aug. 1972 


STREAHTUBE CURVATURE PROGRAM 


Page of 

STC/sheet-2 


Boundary Coordinates 



STREAMTUBE CURVATURE PROGRAM 

Feb. 1972 Page ■ of 

Channel Flow Properties STC/Sheet-3 


— 

channel 



* 

name 




W ' w 




3 CHN 




$A 




ratio of 


flow rate 


specific 

gas 

may be 


heats 

constant 

adjusted? 


(1.4) 

(1.0) 

(T) or F 


GAM= 

, RG= , VARY= 


stagnation properties, see notes 3 and 4 

— 

total temp 

total pressure 



TT0= 

, PT0= 

» 


Mach no. 

static temp 

static pressure 



MACH0= 

. TS0= 

PS0= 


flow area 

normalized by A^ L 



A0= 

» 



$ 


202 


General Instruction and Notes for Sheet-1 of the STC Input Forms 

1) The STC Program computes the subsonic and transonic field of inviscid flow 
past (and within) arbitrarily shaped planar and axisymmetric bodies. Inlet 
and exhaust nozzle flows wherein there may exist jet streams with differing 
energies are typical applications. 

2) The total flow is composed of one or more streams, the properties of 

which are to be listed on Sheet-3 (except as noted below). Each stream 
occupies a "channel" which is Identified by a one to six character 
alphanumeric word. Each channel must be bounded, at least in part, by 
an "upper boundary" and a "lower boundary". Each boundary is also given 
an identifying one to six character name and the coordinates are listed 
on Sheet-2. The following sketch illustrates the naming of channels 
(CHN) and boundaries (BDY) . 



An external flow channel must be named EXT, the recommended name for 
the inlet capture flow channel is INT, and the far-field interface 
boundary must be named FF. Otherwise the selection of the channel and 
boundary names is arbitrary. The special channel names EXT and INT cause 
extra streamlines to be placed in the first refined grid. The boundary 
name FF indicates that the boundary condition on FF is to be obtained 
from an analytic far-field solution. 

There is no specific limit to the complexity of the flow field in regard 
to the number of channels or the number of boundaries. Limits are set 
on the total amount of data which may be input. 
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3) The solution method consists of constructing a grid of streamlines and 
orthogonal lines. Starting with two streamlines per channel (one for 
each boundary) and an orthogonal passing through the first and last 
point of each boundary, the grid is automatically refined by dividing 
the grid intervals in half and in half again as required. The numerical 
resolution, the solution accuracy and the computer execution cost are 
all directly related to the extent of grid refinement. The input variable 
MAXIT determines the maximum number of refinements. Providing this limit 
is not exceeded, the grid will be refined, locally as required, until 
the spacing of orthogonals and streamlines is less than the value deter- 
mined from the SGR and SGZ tables and the Mach number difference between 
any two points on a streamline or an orthogonal line is less than VMG1 
and VMG2, respectively. Grid sire values versus radius (or y-ordinate) 
are to be tabulated after SGR and GR, respectively. NGR is the number 
of entries in each list. Grid sire versus the axial coordinate is to 
be tabulated after SGZ and GZ , respectively. NGZ is the number of GZ 
values. If dimensional values of RG, TSO and PSO are input (see Note 6), 
then VMG1 and VMG2 must have units of velocity rather than Mach number. 

See supplemental notes for additional details. 

A partially refined grid may be saved on tape by specifying a T in 
column 24 of the first card, or read from a previously created tape 
by specifying at T in column 14. 

If TAPE 1 and/or TAPE 2 are not assigned via a REQUEST card, they 
are assigned to disc. This allows the user to obtain output for a 
given refinement level and provides the option of changing input 
parameters on the restart. For the restart case, specify a T in 
column 14 of the first data card and include in the $A list only 
those input quantities (viz; MAXIT) which differ from those originally 
input . 

4. In the initial calculation grid, an orthogonal line will pass through 

each leading and trailing edge point and through each sharp comer point 
(i.e. a point on the boundary with an angle discontinuity). It is not 
possible to analyze a configuration in which two or more of these points 
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are approximately opposite to each other. For example, if a configuration 
contains more than one leading edge, the edges must be staggered relative 
to the streamwise direction. 

5. A free stream Mach number is specified by supplying a value of MACH0. 

6. Perfect gas assumptions are employed and the levels of ambient pressure 
and temperature may be dimensionless (TS0=PS0=RG=1) or dimensional. 
Dimensional values in a consistent set of units as described in 
Reference 5 should be supplied if boundary layer calculations are requested. 

7. A reference (or highlight) area is calculated from the input value of 


R^ as follows: 

axisymmetric : 

= 71 

planar: V ' 4y HL ' V 

This reference area (or Ay in the planar case) is used for defining the 
mass flow for each channel. See STC/Sheet-3 note . 

Computed pressure 
where 

drag forces are normalized by the (maximum) body area 

axisymmetric: 

A „ 2 

Am = tt R 

m 

planar: 

A = Ay = R 
m 'm m 


9. Finite trailing edge thickness is permitted; the maximum thickness, 
or body closure tolerance, is to be supplied after TTE. 

10. On this and the following input sheets, the values in parenthesis are 
used if other values are not input. 

Notes for Sheet-2 of the STC Input Forms 

1. Use one of these sheets for each boundary. Supply a one to six character 
name to identify the boundary in column 14 of the first card. Also 
indicate the name of the channel to which the boundary is adjacent in 
column 24. On the second card indicate whether the boundary is above 
(UPPER = T) or below (UPPER = F) the channel. 
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2. 


The upper or lower "contour" which bounds a given stream may be composed 
of several ’'boundaries". In this case an input sheet must be completed 
for each boundary; the last point of the first boundary must have the 
same coordinates as the first point of the second boundary, and so forth. 

This option is useful when considering variable geometry configurations 
such as flaps or movable nozzle parts. The movable part may be translated 
and rotated, as indicated by Note 8, while the fixed part is held stationary. 

3. List values of Z (or X) , R (or Y) and the surface angle in degrees at 
discrete points along this boundary contour after the symbol "B(l)=". 

Points at sharp comers must be listed twice, one time for each angle 
which exists at that point. In each interval, the STC Program fits 

a locally rotated cubic polynomial. The input points must be smooth 
and consistent with the specified angles. 

All points are to be listed in the streamwise direction. For an inlet 
lip, the points are listed by starting at the highlight point and then 
proceeding around the nose to the trailing edge or downstream boundary. 

The internal and external surfaces are listed separately under different 
boundary names. However, the coordinates of the first point must be the 
same with AN CD equal to +90° for the external surface and -90° for 
the internal surface. 

It is recommended that the boundary coordinates and angles be obtained 
from an analytic definition of the contour, and that around the nose, 
angle variations between points be 20° or less. 

4. Pressure and Mach number distribution data will be printed at each 
orthogonal intersection with the boundary, and not at each input boundary 
point. Orthogonal stations, however, will be placed at any repeated 
point in the boundary table. List the same points twice if it is desired 
to have an orthogonal placed in that position. (This option is modified 
slightly when ZRONLY = T.) Orthogonal stations are always placed at the 
beginning and end of a boundary and at a juncture point between boundaries 
along the same contour. 

5. If the coordinates but not the angles are known, the third column in the 
B-input array may be omitted. In this case specify ZRONLY = T and list the 
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coordinates twice at any point where a curvature jump or an angle jump 
exists. The double points will later be deleted if the angle discontinuity 
is less than 0.01 degrees. These double points are removed because extra 
calculation stations (see Note 5) are usually not desired at such points. 
However, the double point angle tolerance, DBLPTS, preset as 0.01, may be 
input as zero if such double points are to be retained. 


6. With either input option, care must be taken to specify the coordinates 

with precision. The round off or reading error of the coordinate data 

2 

should be less than AS /(10*L), where AS is the local distance between 

points and L is some characteristic length, say the length of the cowl. 

1/2 

Conversely, the spacing between points should not be less than (10 6L) 
where 6 is the relative accuracy of the coordinate data. The tabulated 
output curvatures may be consulted to verify the smoothness of the input 
data. 

7. NACA Series 1 Cowl coordinates are stored internally. With the ZR0NLY = T 
option they may be selected by listing: 

Series 1 Segment • 

Cowl Aft Segment 

where X^, are the highlight coordinates and , Y^ is the position 
of the maximum diameter at the end of the Series 1 contour segment. 


B( 1) = 991, 1, 



8. The input coordinates of a boundary may be adjusted by supplying the 
following input quantities not shown on the front of this sheet: 


ROTATE 

ZPIVOT 

RPIV0T 

SCALE 

ZTRANS 

RTRANS 


angular rotation in degrees 
pivot coordinates 

multiplicative constant to be applied to the coordinate data 

translation increment in the axial direction 

translation on increment in the radial/ vertical direction 
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The order of transformation is rotation, scaling and translation. Hence, 
the pivot coordinates are in the same coordinate frame as the input data 
and the translation increments are in the rotated coordinate frame after 
scaling. It is only necessary to input data for the transformation 
operations to be executed. 

Notes for Sheet-3 of the STC Input Forms 

1. Use one of these sheets for each channel to supply entrance flow properties. 
(See exception under Note 5). 

2. Of the input items shown on the face of the input sheet, use only those 
which are required for the selected options. 

3. The total pressure and total temperature may be input by either of the 
following two procedures: 

a) Specify TTO and PTO if the stagnation properties are known. 

These values may be normalized by the free stream ambient 
temperature and pressure . 

b) Specify MACHO, TSO and PSO if the static properties and Mach 
number are known. Again TSO and PSO may be normalized by the 
free stream ambient values. If only MACHO is supplied (TSO & 

PSO are omitted) the TSO and PSO values from STC/Sheet-1 will 
be used. 

If neither of the above is input, free stream values as supplied on 
Sheet-1 are used for MACHO, PTO and TTO. 

4. If the gas constant, RG, is different from the value supplied on 
STC/Sheet-1, supply the value which applies to this channel. RG, TS0 
TTO, PSO and PTO may all be given as dimensionless (normalized by free 
stream ambient) or dimensional using a consistant set of units. 

5. Input a value AO for the determination of the channel flow rate. AO is 
an area fraction mormalized by A^ as defined under Note 7 of Sheet-1; 
the (dimensional) channel flow area is then the product AO*A HT . The 
flow rate for the channel is computed by using one-dimensional relations 
from the total properties (as determined under Notes 3 & 4) , the 
supplied Mach number, MACH0, and the flow area. For internal inlet 
channels, specify R^ as the highlight radius and AO as the mass flow 
ratio. 
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5 . 


If for any channel the Input data on this sheet is not supplied, the 
reference properties on STC/Sheet-1 will be employed and the frontal 
area calculated at the entrance station will be used. This option is 
suggested for an external stream. 

7. Although approximate flow rates must always be supplied according to 
Note 5, the program will adjust channel flow rates as required 
to meet the zero pressure loading conditions at a trailing edge or to 
meet a maximum (choked) flow rate. The number of channels which require 
flow rate adjustment is equal to the number of trailing edges. If the 
flow rate is not to be varied for this channel, specify VARY = F. 

Grid Refinement Criteria - Supplemental Notes 

One of the important controls that the user should exercise in operating 
the STC program is the Grid Refinement Specification. In simpliest terms the 
user is required to input the following three quantities: 

SGR = , Nominal distance between streamlines and between 

orthogonal lines 

VMG1 = - , Nominal Mach number (velocity) difference between 

points along a streamline 

VMG2 = , Nominal Mach number (velocity) difference between 

points along an orthogonal line 

Note, if dimensional values of the perfect gas constant and temperature are 
input then VMG1 and VMG2 are velocity differences; otherwise VMG1 and VMG2 
are to be input as a Mach number difference. 

To illustrate the function of these input items, consider the question of 
whether or not a new orthogonal line should be inserted between existing 
orthogonals as shown in Figure la. To answer this question, the program 
computes, say, between the (i) and (i+1) orthogonals, the following quantities 


at each 

streamline: 



(la) 

AS1 

R = 

SG : 

(SG - SGR) 

(la) 
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(2a) 



1 - V 

JL+Jhi Lti 

VMGl 


j - 1, 2, 3,.. . (lb) 


Figuratively speaking, these values are then plotted on Figure 2, 


If all the RSj ,Rfj points fall outside of Region (1) a new orthogonal will not 
be inserted into the i,i+l interval. But if one or more points do fall in 
Region (1), a new orthogonal is required. The new line will be extended to 
all streamlines for which the rs^.rVj values fall in Region (2), or to a 
minimum of five streamlines. 


Thus, according to the criteria illustrated in Figure 2, the maximum distance 
between any two points along a streamline is SG (=SGR) . If a velocity gradient 
exists along the streamline, the maximum distance between points is reduced. 
When the spacing ratio Rs becomes very small (Rs < 0.24) relatively large 
velocity differences between points is allowed, a feature which is intended 
to eliminate excessive refinement near stagnation points. 


A similar criteria is used for determining when additional streamlines are 
required. The only variation is that Eqs. (la) and (2a) are replaced by 
Eqs . (lb) and (2b): 


v. 

l 


AS 2 
SG 

|v i-m - Yj 

VMG2 


1, 2, 3... 


(lb) 

(2b) 


Streamlines and orthogonal lines are inserted until the above criteria is 
satisfied (or until the number of refinements equals MAXIT) . The internal 
values of VMGl and VMG2 (0.1) are appropriate for most configurations, and 
unless these values are to be altered the VMGl and VMG2 input may be omitted. 
However, the user must always supply a value for SGR. In selecting this value 
care must be taken to ensure that a moderate grid will develop so that reasonable 
accuracy is obtained without sacrificing computing efficiency with an excessive 
number of grid points. * 

As a guide to the selection of SGR it is suggested that SGR be set to about 
twice the expected grid spacing. Because new grid intervals are always obtained 
by halving existing grid intervals, the average grid size will be less than 
that required to satisfy the criteria given in Figure 2. 
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Spatial Variation of Grid Size 


For external flow problems, and some internal problems as well, it is necessary 
to input a spatial grid size variation. The additional input required is: 

GR = » » > > » » 

SGR - > t , » > > 

= > > » * » > 

SGZ = , • • . % . 


GR and SGR is a table of grid size (SGR) with radius (GR) ; likewise GZ and 
SGZ is a table of grid size (GZ) with axial position (SGZ). The value of SG 
in Eqs . (la) and (lb) at any point, z,r, is then taken as the maximum of the 
SG's found by linear interpolation in the two tables of radial and axial 
variations. 


Special rules which apply to these input items are as follows: 

1) If a radial variation (but no axial variation) of grid size is desired, 
the user must input pairs of values in the SGR and GR lists. GR must 
be in ascending order. If grid size information is required beyond the 
limits of the input table, the closest values will be used. Also if 
the grid size is to be constant with radius no GR input is required and 
only one value of SGR need be supplied. 

2) If an axial variation (but no radial variation) of grid size is desired, 
the user must supply pairs of values in the SGZ and GZ lists. Also he 
must input one value of SGR which is thy minimum of the SGZ values. 

3) If both an axial and radial variation of grid size is desired, he should 
supply appropriate values in both tables. Remember that for any point, 
r,z, the r-coordinate will be used to obtain a grid size from the GR.SGR 
table and the z-coordinate will be used to interpolate in the GZ,SGZ table 
The maximum of the two interpolated values will then be used in Eqs. (la) 
and (lb). Generally, the SGR and SGZ values are set to small values close 
to a body, or region of interest, and increased to larger values as one 
moves away from the body. 
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Input Control of "Partial" Grid Lines 1 

As noted previously, a new orthogonal lines and streamlines are not necessarily 
extended to the field boundary. Instead, the length of these lines is determined 
by the sire of Region (2) in Figure 2. Notice that Region (2) is increased, and 
the length of the new grid lines is likewise increased, if CRX, the constant in 
the equation defining the lower boundary of the region, is decreased. If 
CRX = 0, then all new grid lines will span the field; that is, full grid lines 
rather than "partial" grid lines will be inserted. 

Values of CRX are preset, but they may be altered by input cards as follows: 


Input Name 

Preset 

Used for Extending 

(a) 

(b) 

Value 


CRX(l) 

CRXSL 

.375 

New Streamlines 

CRX (2) 

CRXOL 

.375 

New Orthogonal lines across 
a subsonic region 

CRX(3) 

CRXSS 

.125 

New Orthogonal lines across 
a supersonic or mixed flow 
region 

CRX(4) 

CRXE 

0 

New orthogonal lines which 
cross a sonic line 

CRX(5) 

CRXC 

0 

New Orthogonal lines which 
cross a supersonic to subsonic 
compression line 

two or more of 

the above values 

app ly , the 

smaller value will be used. 

Recommended 

Spatial Grid Variation for an 

Inlet Configuration 


Typical input for the spatial variation of grid size for an inlet configuration 
is given in Figure 3. In this case all dimensions are normalized by the highlight 
radius. The quoted values of SGR, GR, SGZ, GZ should be scaled by the highlight 
radius for a given problem. Notice that the combination of the axial and radial 
variations yields contours of grid size which are roughly centered about the 
leading edge. 
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Fig. 12-2 Grid Refinement Criteria 
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